• Monday: First lecture: Presentation of the course, aims and content
• Monday: Second Lecture: Introduction to C++ programming and numerical precision.
• Tuesday: Numerical precision and C++ programming, continued and exercises for first week
• Numerical differentiation and loss of numerical precision (chapter 3 lecture notes)
• Computer lab: Thursday and Friday. First time: Thursday and Friday this week, Presentation of hardware and software at room FV329 first hour of every labgroup and solution of first simple exercises.

• Lectures: Monday (4.15pm-6pm) and Tuesday (12.15pm-2pm) only this and next week. Thereafter weeks 38 and 39, and weeks 43 and 44 and finally weeks 47 and 48.
• Weekly reading assignments needed to solve projects.
• First hour of each lab session used to discuss technicalities, address questions etc linked with projects.
• Detailed lecture notes, exercises, all programs presented, projects etc can be found at the homepage of the course.
• Computerlab: Thursday (9am-7pm) and Friday (9am-7pm) room FV329.
• Weekly plans and all other information are on the official webpage.
• Final written exam December 12, 9am (four hours).

• Several computer exercises, 5 compulsory projects. Electronic reports only.
• Evaluation and grading: The last project (50% of final grade) and a final written exam (50% of final grade). Final written exam December 12.
• The computer lab (room FV329)consists of 16 Linux PCs, but many prefer own laptops. C/C++ is the default programming language, but Fortran2008 and Python are also used. All source codes discussed during the lectures can be found at the webpage of the course. We recommend either C/C++, Fortran2008 or Python as languages.

day	teacher
Thursday 9am-1pm	Anders, Morten L., Håvard, MHJ
Thursday 1pm-5pm	Anders, Morten L., Håvard, MHJ
Friday 9am-1pm	Anders, Morten L., Håvard, MHJ
Friday 1pm-5pm	Anders, Morten L., Håvard, MHJ

• Numerical precision and intro to C++ programming
• Numerical derivation and integration
• Random numbers and Monte Carlo integration
• Monte Carlo methods in statistical physics
• Quantum Monte Carlo methods
• Linear algebra and eigenvalue problems
• Non-linear equations and roots of polynomials
• Ordinary differential equations
• Partial differential equations
• Parallelization of codes
• Programming av GPUs (optional)

Linear algebra and eigenvalue problems, chapters 6 and 7.

• Know Gaussian elimination and LU decomposition
• How to solve linear equations
• How to obtain the inverse and the determinant of a real symmetric matrix
• Cholesky and tridiagonal matrix decomposition

Linear algebra and eigenvalue problems, chapters 6 and 7.

• Householder's tridiagonalization technique and finding eigenvalues based on this
• Jacobi's method for finding eigenvalues
• Singular value decomposition
• Qubic Spline interpolation

Numerical integration, standard methods and Monte Carlo methods (chapters 4 and 11).

• Trapezoidal, rectangle and Simpson's rules
• Gaussian quadrature, emphasis on Legendre polynomials, but you need to know about other polynomials as well.
• Brute force Monte Carlo integration
• Random numbers (simplest algo, ran0) and probability distribution functions, expectation values
• Improved Monte Carlo integration and importance sampling.

Monte Carlo methods in physics (chapters 12, 13, and 14).

• Random walks and Markov chains and relation with diffusion equation
• Metropolis algorithm, detailed balance and ergodicity
• Simple spin systems and phase transitions
• Variational Monte Carlo
• How to construct trial wave functions for quantum systems

Ordinary differential equations (chapters 8 and 9).

• Euler's method and improved Euler's method, truncation errors
• Runge Kutta methods, 2nd and 4th order, truncation errors
• How to implement a second-order differential equation, both linear and non-linear. How to make your equations dimensionless.
• Boundary value problems, shooting and matching method (chap 9).

Partial differential equations, chapter 10.

• Set up diffusion, Poisson and wave equations up to 2 spatial dimensions and time
• Set up the mathematical model and algorithms for these equations, with boundary and initial conditions. Their stability conditions.
• Explicit, implicit and Crank-Nicolson schemes, and how to solve them. Remember that they result in triangular matrices.
• How to compute the Laplacian in Poisson's equation.
• How to solve the wave equation in one and two dimensions.

Bachelor/Master/PhD Courses.

• INF-MAT4350 Numerical linear algebra
• MAT-INF3300/3310, PDEs and Sobolev spaces I and II
• INF-MAT3360 PDEs
• INF5620 Numerical methods for PDEs, finite element method
• FYS4411 Computational physics II (Parallelization (MPI), object orientation, quantum mechanical systems with many interacting particles), spring semester
• FYS4460 Computational physics III (Parallelization (MPI), object orientation, classical statistical physics, simulation of phase transitions, spring semester
• INF3331 Problem solving with high-level languages (Python), fall semester
• INF3380 Parallel computing for problems in the Natural Sciences (mostly PDEs), spring semester

and discussed at the lab sessions the first week.

• GIT for version control (see webpage)
• ipython notebook
• QTcreator for editing and mastering computational projects (for C++ codes, see webpage of course)
• Armadillo as a useful numerical library for C++, highly recommended
• Unit tests, see also webpage
• Devilry for handing in projects

Compiling and linking, without QTcreator.

In order to obtain an executable file for a C++ program, the following instructions under Linux/Unix can be used

c++ -c -Wall myprogram.cpp
c++ -o myprogram myprogram.o

where the compiler is called through the command c++/g++. The compiler option -Wall means that a warning is issued in case of non-standard language. The executable file is in this case myprogram. The option -c is for compilation only, where the program is translated into machine code, while the -o option links the produced object file myprogram.o and produces the executable myprogram .

For Fortran2008 we use the Intel compiler, replace c++ with ifort. Also, to speed up the code use compile options like

c++ -O3 -c -Wall myprogram.cpp

# Comment lines
# General makefile for c - choose PROG =   name of given program
# Here we define compiler option, libraries and the  target
CC= g++ -Wall
PROG= myprogram
# this is the math library in C, not necessary for C++
LIB = -lm
# Here we make the executable file
${PROG} :          ${PROG}.o
                   ${CC} ${PROG}.o ${LIB} -o ${PROG}
# whereas here we create the object file
${PROG}.o :       ${PROG}.c
                  ${CC} -c ${PROG}.c

The C encounter.

Here we present first the C version.

/* comments in C begin like this and end with */
#include <stdlib.h> /* atof function */
#include <math.h>   /* sine function */
#include <stdio.h>  /* printf function */
int main (int argc, char* argv[])
{
  double r, s;        /* declare variables */
  r = atof(argv[1]);  /* convert the text argv[1] to double */
  s = sin(r);
  printf("Hello, World! sin(%g)=%g\n", r, s);
  return 0;           /* success execution of the program */

Dissection I.

The compiler must see a declaration of a function before you can call it (the compiler checks the argument and return types). The declaration of library functions appears in so-called "header files" that must be included in the program, e.g.,

   #include <stdlib.h> /* atof function */

We call three functions (atof, sin, printf) and these are declared in three different header files. The main program is a function called main with a return value set to an integer, int (0 if success). The operating system stores the return value, and other programs/utilities can check whether the execution was successful or not. The command-line arguments are transferred to the main function through

   int main (int argc, char* argv[])

Dissection II.

The command-line arguments are transferred to the main function through

   int main (int argc, char* argv[])

The integer argc is the no of command-line arguments, set to one in our case, while argv is a vector of strings containing the command-line arguments with argv[0] containing the name of the program and argv[1], argv[2], ... are the command-line args, i.e., the number of lines of input to the program. Here we define floating points, see also below, through the keywords float for single precision real numbers and double for double precision. The function atof transforms a text (argv[1]) to a float. The sine function is declared in math.h, a library which is not automatically included and needs to be linked when computing an executable file.

With the command printf we obtain a formatted printout. The printf syntax is used for formatting output in many C-inspired languages (Perl, Python, Awk, partly C++).

Now in C++.

Here we present first the C++ version.

// A comment line begins like this in C++ programs
// Standard ANSI-C++ include files
using namespace std
#include <iostream>  // input and output
int main (int argc, char* argv[])
{
  // convert the text argv[1] to double using atof:
  double r = atof(argv[1]);
  double s = sin(r);
  cout << "Hello, World! sin(" << r << ")=" << s << '\n';
  // success
  return 0;
}

Dissection I.

We have replaced the call to printf with the standard C++ function cout. The header file <iostream.h> is then needed. In addition, we don't need to declare variables like r and s at the beginning of the program. I personally prefer however to declare all variables at the beginning of a function, as this gives me a feeling of greater readability.

C/C++ program.

• A C/C++ program begins with include statements of header files (libraries,intrinsic functions etc)
• Functions which are used are normally defined at top (details next week)
• The main program is set up as an integer, it returns 0 (everything correct) or 1 (something went wrong)
• Standard if, while and for statements as in Java, Fortran, Python...
• Integers have a very limited range.

Arrays.

• A C/C++ array begins by indexing at 0!
• Array allocations are done by size, not by the final index value.If you allocate an array with 10 elements, you should index them from \( 0,1,\dots, 9 \).
• Initialize always an array before a computation.

Integer and Real Numbers.

• Overflow
• Underflow
• Roundoff errors
• Loss of precision

C++ and Fortran declarations.

type in C/C++ and Fortran2008	bits	range
int/INTEGER (2)	16	-32768 to 32767
unsigned int	16	0 to 65535
signed int	16	-32768 to 32767
short int	16	-32768 to 32767
unsigned short int	16	0 to 65535
signed short int	16	\( -32768 \) to 32767
int/long int/INTEGER (4)	32	-2147483648 to 2147483647
signed long int	32	-2147483648 to 2147483647
float/REAL(4)	32	\( 3.4\times 10^{-44} \) to \( 3.4\times 10^{+38} \)
double/REAL(8)	64	\( 1.7\times 10^{-322} \) to \( 1.7\times 10^{+308} \)
long double	64	\( 1.7\times 10^{-322} \) to \( 1.7\times 10^{+308} \)

How to do it.

$$ a_n2^n+a_{n-1}2^{n-1} +a_{n-2}2^{n-2} +\dots +a_{0}2^{0}. $$

In binary notation we have thus \( (417)_{10} =(110110001)_2 \) since we have $$ (110100001)_2 =1\times2^8+1\times 2^{7}+0\times 2^{6}+1\times 2^{5}+0\times 2^{4}+0\times 2^{3}+0\times 2^{2}+0\times 2^{2}+0\times 2^{1}+1\times 2^{0}. $$ To see this, we have performed the following divisions by 2

417/2=208	remainder 1	coefficient of \( 2^{0} \) is 1
208/2=104	remainder 0	coefficient of \( 2^{1} \) is 0
104/2=52	remainder 0	coefficient of \( 2^{2} \) is 0
52/2=26	remainder 0	coefficient of \( 2^{3} \) is 0
26/2=13	remainder 1	coefficient of \( 2^{4} \) is 0
13/2= 6	remainder 1	coefficient of \( 2^{5} \) is 1
6/2= 3	remainder 0	coefficient of \( 2^{6} \) is 0
3/2= 1	remainder 1	coefficient of \( 2^{7} \) is 1
1/2= 0	remainder 1	coefficient of \( 2^{8} \) is 1

Integer numbers.

using namespace std;
#include <iostream>
int main (int argc, char* argv[])
{
  int i;
  int terms[32]; // storage of a0, a1, etc, up to 32 bits
  int number = atoi(argv[1]);
  // initialise the term a0, a1 etc
  for (i=0; i < 32 ; i++){ terms[i] = 0;}
  for (i=0; i < 32 ; i++){
    terms[i] = number%2;
    number /= 2;

  // write out results
  cout << "Number of bytes used= " << sizeof(number) << endl;
  for (i=0; i < 32 ; i++){
    cout << " Term nr: " << i << "Value= " << terms[i];
    cout << endl;

  return 0;

Integer numbers, Fortran.

PROGRAM binary_integer
IMPLICIT NONE
  INTEGER  i, number, terms(0:31) ! storage of a0, a1, etc, up to 32 bits

  WRITE(*,*) 'Give a number to transform to binary notation'
  READ(*,*) number
! Initialise the terms a0, a1 etc
  terms = 0
! Fortran takes only integer loop variables
  DO i=0, 31
     terms(i) = MOD(number,2)
     number = number/2
  ENDDO
! write out results
  WRITE(*,*) 'Binary representation '
  DO i=0, 31
    WRITE(*,*)' Term nr and value', i, terms(i)
  ENDDO

END PROGRAM binary_integer

Possible Overflow for Integers.

// A comment line begins like this in C++ programs
// Program to calculate 2**n
// Standard ANSI-C++ include files */
using namespace std
#include <iostream>
#include <cmath>
int main()
{
   int  int1, int2, int3;
// print to screen
   cout << "Read in the exponential N for 2^N =\n";
// read from screen
   cin >> int2;
   int1 = (int) pow(2., (double) int2);
   cout << " 2^N * 2^N = " << int1*int1 << "\n";
   int3 = int1 - 1;
   cout << " 2^N*(2^N - 1) = " << int1 * int3  << "\n";
   cout << " 2^N- 1 = " << int3  << "\n";
   return 0;

// End: program main()

Machine Numbers.

In the decimal system we would write a number like \( 9.90625 \) in what is called the normalized scientific notation. $$ 9.90625=0.990625\times 10^{1}, $$ and a real non-zero number could be generalized as $$ \begin{equation} x=\pm r\times 10^{{\mbox{n}}}, \end{equation} $$ with \( r \) a number in the range \( 1/10 \le r < 1 \). In a similar way we can use represent a binary number in scientific notation as $$ \begin{equation} x=\pm q\times 2^{{\mbox{m}}}, \end{equation} $$ with \( q \) a number in the range \( 1/2 \le q < 1 \). This means that the mantissa of a binary number would be represented by the general formula $$ \begin{equation} (0.a_{-1}a_{-2}\dots a_{-n})_2=a_{-1}\times 2^{-1} +a_{-2}\times 2^{-2}+\dots+a_{-n}\times 2^{-n}. \end{equation} $$

Machine Numbers.

In a typical computer, floating-point numbers are represented in the way described above, but with certain restrictions on \( q \) and \( m \) imposed by the available word length. In the machine, our number \( x \) is represented as $$ \begin{equation} x=(-1)^s\times {\mbox{mantissa}}\times 2^{{\mbox{exponent}}}, \end{equation} $$

where \( s \) is the sign bit, and the exponent gives the available range. With a single-precision word, 32 bits, 8 bits would typically be reserved for the exponent, 1 bit for the sign and 23 for the mantissa.

Machine Numbers.

A modification of the scientific notation for binary numbers is to require that the leading binary digit 1 appears to the left of the binary point. In this case the representation of the mantissa \( q \) would be \( (1.f)_2 \) and $ 1 \le q < 2$. This form is rather useful when storing binary numbers in a computer word, since we can always assume that the leading bit 1 is there. One bit of space can then be saved meaning that a 23 bits mantissa has actually 24 bits. This means explicitely that a binary number with 23 bits for the mantissa reads $$ \begin{equation} (1.a_{-1}a_{-2}\dots a_{-23})_2=1\times 2^0+a_{-1}\times 2^{-1}+ +a_{-2}\times 2^{-2}+\dots+a_{-23}\times 2^{-23}. \end{equation} $$ As an example, consider the 32 bits binary number $$ (10111110111101000000000000000000)_2, $$ where the first bit is reserved for the sign, 1 in this case yielding a negative sign. The exponent \( m \) is given by the next 8 binary numbers \( 01111101 \) resulting in 125 in the decimal system.

Machine Numbers.

However, since the exponent has eight bits, this means it has \( 2^8-1=255 \) possible numbers in the interval \( -128 \le m \le 127 \), our final exponent is \( 125-127=-2 \) resulting in \( 2^{-2} \). Inserting the sign and the mantissa yields the final number in the decimal representation as $$ -2^{-2}\left(1\times 2^0+1\times 2^{-1}+ 1\times 2^{-2}+1\times 2^{-3}+0\times 2^{-4}+1\times 2^{-5}\right)=$$ $$ (-0.4765625)_{10}. $$ In this case we have an exact machine representation with 32 bits (actually, we need less than 23 bits for the mantissa).

If our number \( x \) can be exactly represented in the machine, we call \( x \) a machine number. Unfortunately, most numbers cannot and are thereby only approximated in the machine. When such a number occurs as the result of reading some input data or of a computation, an inevitable error will arise in representing it as accurately as possible by a machine number.

Machine Numbers.

A floating number x, labelled \( fl(x) \) will therefore always be represented as $$ \begin{equation} fl(x) = x(1\pm \epsilon_x), \end{equation} $$ with \( x \) the exact number and the error \( |\epsilon_x| \le |\epsilon_M| \), where \( \epsilon_M \) is the precision assigned. A number like \( 1/10 \) has no exact binary representation with single or double precision. Since the mantissa $$ \left(1.a_{-1}a_{-2}\dots a_{-n}\right)_2 $$ is always truncated at some stage \( n \) due to its limited number of bits, there is only a limited number of real binary numbers. The spacing between every real binary number is given by the chosen machine precision. For a 32 bit words this number is approximately $ \epsilon_M \sim 10^{-7}$ and for double precision (64 bits) we have $ \epsilon_M \sim 10^{-16}$, or in terms of a binary base as \( 2^{-23} \) and \( 2^{-52} \) for single and double precision, respectively.

Machine Numbers.

In the machine a number is represented as $$ \begin{equation} fl(x)= x(1+\epsilon) \end{equation} $$

where \( |\epsilon| \leq \epsilon_M \) and \( \epsilon \) is given by the specified precision, \( 10^{-7} \) for single and \( 10^{-16} \) for double precision, respectively. \( \epsilon_M \) is the given precision. In case of a subtraction \( a=b-c \), we have $$ \begin{equation} fl(a)=fl(b)-fl(c) = a(1+\epsilon_a), \end{equation} $$ or $$ \begin{equation} fl(a)=b(1+\epsilon_b)-c(1+\epsilon_c), \end{equation} $$

meaning that $$ \begin{equation} fl(a)/a=1+\epsilon_b\frac{b}{a}- \epsilon_c\frac{c}{a}, \end{equation} $$

and if \( b\approx c \) we see that there is a potential for an increased error in \( fl(a) \).

Machine Numbers.

Define the absolute error as $$ \begin{equation} |fl(a)-a|, \end{equation} $$ whereas the relative error is $$ \begin{equation} \frac{ |fl(a)-a|}{a} \le \epsilon_a. \end{equation} $$ The above subraction is thus $$ \begin{equation} \frac{ |fl(a)-a|}{a}=\frac{ |fl(b)-fl(c)-(b-c)|}{a}, \end{equation} $$ yielding $$ \begin{equation} \frac{ |fl(a)-a|}{a}=\frac{ |b\epsilon_b- c\epsilon_c|}{a}. \end{equation} $$ The relative error is the quantity of interest in scientific work. Information about the absolute error is normally of little use in the absence of the magnitude of the quantity being measured.

Real Numbers.

• Overflow: When the positive exponent exceeds the max value, e.g., 308 for DOUBLE PRECISION (64 bits). Under such circumstances the program will terminate and some compilers may give you the warning OVERFLOW.
• Underflow: When the negative exponent becomes smaller than the min value, e.g., -308 for DOUBLE PRECISION. Normally, the variable is then set to zero and the program continues. Other compilers (or compiler options) may warn you with the UNDERFLOW message and the program terminates.

Real Numbers.

• Loss of precision: When one has to e.g., multiply two large numbers where one suspects that the outcome may be beyond the bonds imposed by the variable declaration, one could represent the numbers by logarithms, or rewrite the equations to be solved in terms of dimensionless variables. When dealing with problems in e.g., particle physics or nuclear physics where distance is measured in fm (\( 10^{-15} \) m), it can be quite convenient to redefine the variables for distance in terms of a dimensionless variable of the order of unity. To give an example, suppose you work with single precision and wish to perform the addition \( 1+10^{-8} \). In this case, the information containing in \( 10^{-8} \) is simply lost in the addition. Typically, when performing the addition, the computer equates first the exponents of the two numbers to be added. For \( 10^{-8} \) this has however catastrophic consequences since in order to obtain an exponent equal to \( 10^0 \), bits in the mantissa are shifted to the right. At the end, all bits in the mantissa are zeros.

Three ways of computing \( e^{-x} \).

Brute force: $$\exp{(-x)}=\sum_{n=0}^{\infty}(-1)^n\frac{x^n}{n!}$$

Recursion relation for $$ \exp{(-x)}=\sum_{n=0}^{\infty}s_n=\sum_{n=0}^{\infty}(-1)^n\frac{x^n}{n!} $$ $$ s_n=-s_{n-1}\frac{x}{n}, $$ $$ \exp{(x)}=\sum_{n=0}^{\infty}s_n $$ $$ \exp{(-x)}=\frac{1}{\exp{(x)}} $$

Brute Force.

// Program to calculate function exp(-x)
// using straightforward summation with differing  precision
using namespace std
#include <iostream>
#include <cmath>
// type float:  32 bits precision
// type double: 64 bits precision
#define   TYPE          double
#define   PHASE(a)      (1 - 2 * (abs(a) % 2))
#define   TRUNCATION    1.0E-10
// function declaration
TYPE factorial(int);

Still Brute Force.

int main()
{
   int   n;
   TYPE  x, term, sum;
   for(x = 0.0; x < 100.0; x += 10.0)  {
     sum  = 0.0;                //initialization
     n    = 0;
     term = 1;
     while(fabs(term) > TRUNCATION)  {
         term =  PHASE(n) * (TYPE) pow((TYPE) x,(TYPE) n)
                / factorial(n);
         sum += term;
         n++;
     }  // end of while() loop

Oh, it never ends!

      printf("\nx = %4.1f   exp = %12.5E  series = %12.5E
              number of terms = %d",
              x, exp(-x), sum, n);
   } // end of for() loop

   printf("\n");           // a final line shift on output
   return 0;
} // End: function main()
//     The function factorial()
//     calculates and returns n!
TYPE factorial(int n)
{
   int  loop;
   TYPE fac;
   for(loop = 1, fac = 1.0; loop <= n; loop++)  {
      fac *= loop;

   return fac;
} // End: function factorial()

What is going on?

\( x \)	\( \exp{(-x)} \)	Series	Number of terms in series
0.0	0.100000E+01	0.100000E+01	1
10.0	0.453999E-04	0.453999E-04	44
20.0	0.206115E-08	0.487460E-08	72
30.0	0.935762E-13	-0.342134E-04	100
40.0	0.424835E-17	-0.221033E+01	127
50.0	0.192875E-21	-0.833851E+05	155
60.0	0.875651E-26	-0.850381E+09	171
70.0	0.397545E-30	NaN	171
80.0	0.180485E-34	NaN	171
90.0	0.819401E-39	NaN	171
100.0	0.372008E-43	NaN	171

// program to compute exp(-x) without exponentials
using namespace std
#include <iostream>
#include <cmath>
#define  TRUNCATION     1.0E-10

int main()
{
   int       loop, n;
   double    x, term, sum;
   for(loop = 0; loop <= 100; loop += 10)
   {
     x    = (double) loop;          // initialization
     sum  = 1.0;
     term = 1;
     n    = 1;

Last statements.

     while(fabs(term) > TRUNCATION)
       {
	 term *= -x/((double) n);
	 sum  += term;
	 n++;
       } // end while loop
     cout << "x = " << x   << " exp = " << exp(-x) <<"series = "
          << sum  << " number of terms =" << n << "\n";
   } // end of for() loop

   cout << "\n";           // a final line shift on output

}  /*    End: function main() */

More Problems.

\( x \)	\( \exp{(-x)} \)	Series	Number of terms in series
0.000000	0.10000000E+01	0.10000000E+01	1
10.000000	0.45399900E-04	0.45399900E-04	44
20.000000	0.20611536E-08	0.56385075E-08	72
30.000000	0.93576230E-13	-0.30668111E-04	100
40.000000	0.42483543E-17	-0.31657319E+01	127
50.000000	0.19287498E-21	0.11072933E+05	155
60.000000	0.87565108E-26	-0.33516811E+09	182
70.000000	0.39754497E-30	-0.32979605E+14	209
80.000000	0.18048514E-34	0.91805682E+17	237
90.000000	0.81940126E-39	-0.50516254E+22	264
100.000000	0.37200760E-43	-0.29137556E+26	291

3 point formulae.

First derivative (\( f_0 = f(x_0) \), \( f_{-h}=f(x_0-h) \) and \( f_{+h}=f(x_0+h) \) $$ \frac{f_h-f_{-h}}{2h}=f'_0+\sum_{j=1}^{\infty}\frac{f_0^{(2j+1)}}{(2j+1)!}h^{2j}. $$ Second derivative $$ \frac{ f_h -2f_0 +f_{-h}}{h^2}=f_0''+2\sum_{j=1}^{\infty}\frac{f_0^{(2j+2)}}{(2j+2)!}h^{2j}. $$

$$ \epsilon=log_{10}\left(\left|\frac{f''_{\mbox{computed}}-f''_{\mbox{exact}}} {f''_{\mbox{exact}}}\right|\right), $$ $$ \epsilon_{\mbox{tot}}=\epsilon_{\mbox{approx}}+\epsilon_{\mbox{ro}}. $$

For the computed second derivative we have $$ f_0''=\frac{ f_h -2f_0 +f_{-h}}{h^2}-2\sum_{j=1}^{\infty}\frac{f_0^{(2j+2)}}{(2j+2)!}h^{2j}, $$ and the truncation or approximation error goes like $$ \epsilon_{\mbox{approx}}\approx \frac{f_0^{(4)}}{12}h^{2}. $$

If we were not to worry about loss of precision, we could in principle make \( h \) as small as possible. However, due to the computed expression in the above program example $$ f_0''=\frac{ f_h -2f_0 +f_{-h}}{h^2}=\frac{ (f_h -f_0) +(f_{-h}-f_0)}{h^2}, $$ we reach fairly quickly a limit for where loss of precision due to the subtraction of two nearly equal numbers becomes crucial.

If \( (f_{\pm h} -f_0) \) are very close, we have \( (f_{\pm h} -f_0)\approx \epsilon_M \), where \( |\epsilon_M|\le 10^{-7} \) for single and \( |\epsilon_M|\le 10^{-15} \) for double precision, respectively.

We have then $$ \left|f_0''\right|= \left|\frac{ (f_h -f_0) +(f_{-h}-f_0)}{h^2}\right|\le \frac{ 2 \epsilon_M}{h^2}. $$

Our total error becomes $$ \left|\epsilon_{\mbox{tot}}\right|\le \frac{2 \epsilon_M}{h^2} + \frac{f_0^{(4)}}{12}h^{2}. \label{eq:experror} $$ It is then natural to ask which value of \( h \) yields the smallest total error. Taking the derivative of \( \left|\epsilon_{\mbox{tot}}\right| \) with respect to \( h \) results in $$ h= \left(\frac{ 24\epsilon_M}{f_0^{(4)}}\right)^{1/4}. $$ With double precision and \( x=10 \) we obtain $$ h\approx 10^{-4}. $$ Beyond this value, it is essentially the loss of numerical precision which takes over.

Due to the subtractive cancellation in the expression for \( f'' \) there is a pronounced detoriation in accuracy as \( h \) is made smaller and smaller.

It is instructive in this analysis to rewrite the numerator of the computed derivative as $$ (f_h -f_0) +(f_{-h}-f_0)=(e^{x+h}-e^{x}) + (e^{x-h}-e^{x}), $$ as $$ (f_h -f_0) +(f_{-h}-f_0)=e^x(e^{h}+e^{-h}-2), $$ since it is the difference \( (e^{h}+e^{-h}-2) \) which causes the loss of precision.

\( x \)	\( h=0.01 \)	\( h=0.001 \)	\( h=0.0001 \)	\( h=0.0000001 \)	Exact
0.0	1.000008	1.000000	1.000000	1.010303	1.000000
1.0	2.718304	2.718282	2.718282	2.753353	2.718282
2.0	7.389118	7.389057	7.389056	7.283063	7.389056
3.0	20.085704	20.085539	20.085537	20.250467	20.085537
4.0	54.598605	54.598155	54.598151	54.711789	54.598150
5.0	148.414396	148.413172	148.413161	150.635056	148.413159

The results for \( x=10 \) are shown in the Table

\( h \)	\( e^{h}+e^{-h} \)	\( e^{h}+e^{-h}-2 \)
\( 10^{-1} \)	2.0100083361116070	\( 1.0008336111607230\times 10^{-2} \)
\( 10^{-2} \)	2.0001000008333358	\( 1.0000083333605581\times 10^{-4} \)
\( 10^{-3} \)	2.0000010000000836	\( 1.0000000834065048\times 10^{-6} \)
\( 10^{-5} \)	2.0000000099999999	\( 1.0000000050247593\times 10^{-8} \)
\( 10^{-5} \)	2.0000000001000000	\( 9.9999897251734637\times 10^{-11} \)
\( 10^{-6} \)	2.0000000000010001	\( 9.9997787827987850\times 10^{-13} \)
\( 10^{-7} \)	2.0000000000000098	\( 9.9920072216264089\times 10^{-15} \)
\( 10^{-8} \)	2.0000000000000000	\( 0.0000000000000000\times 10^{0} \)
\( 10^{-9} \)	2.0000000000000000	\( 1.1102230246251565\times 10^{-16} \)
\( 10^{-10} \)	2.0000000000000000	\( 0.0000000000000000\times 10^{0} \)

• Monday: Repetition from last week
• Numerical differentiation
• C/C++ programming details, pointers, read/write to/from file
• Tuesday: Intro to linear Algebra and presentation of project 1.
• Matrices in C++ and Fortran2008
• Gaussian elimination and discussion of project 1.
• Computer-Lab: start project 1. Reading asssignments and preparation for project 1: sections 2.5 and 3.1 for general C++ and Fortran features. Sections 6.1-6.4 (till page 182) are relevant for project 1.

A pointer specifies where a value resides in the computer's memory (like a house number specifies where a particular family resides on a street).

A pointer points to an address not to a data container of any kind!

Simple example declarations:

  using namespace std; // note use of namespace
  int main()
 {
   // what are the differences?
   int var;
   cin >> var;
   int *p, q;
   int *s, *t;
   int * a new[var];    // dynamic memory allocation
   delete [] a;
}

using namespace std; // note use of namespace
int main()
{
  int var;
  int *p;
  p = &var;
  var  = 421;
  printf("Address of integer variable var : %p\n",&var);
  printf("Its value: %d\n", var);
  printf("Value of integer pointer p : %p\n",p);
  printf("The value p points at :  %d\n",*p);
  printf("Address of the pointer p : %p\n",&p);
  return 0;
}

Discussion.

int main()
{
  int var;     // Define an integer variable var
  int *p;      // Define a pointer to an integer
  p = &var;    // Extract the address of var
  var = 421;   // Change content of var
  printf("Address of integer variable var : %p\n", &var);
  printf("Its value: %d\n", var);  // 421
  printf("Value of integer pointer p : %p\n", p);  // = &var
  // The content of the variable pointed to by p is *p
  printf("The value p points at :  %d\n", *p);
  // Address where the pointer is stored in memory
  printf("Address of the pointer p : %p\n", &p);
  return 0;
}

int matr[2];
int *p;
p = &matr[0];
matr[0] = 321;
matr[1] = 322;
printf("\nAddress of matrix element matr[1]: %p",&matr[0]);
printf("\nValue of the  matrix element  matr[1]; %d",matr[0]);
printf("\nAddress of matrix element matr[2]: %p",&matr[1]);
printf("\nValue of the matrix element  matr[2]: %d\n", matr[1]);
printf("\nValue of the pointer p: %p",p);
printf("\nThe value p points to: %d",*p);
printf("\nThe value that (p+1) points to  %d\n",*(p+1));
printf("\nAddress of pointer p : %p\n",&p);

int matr[2];    // Define integer array with two elements
int *p;         // Define pointer to integer
p = &matr[0];   // Point to the address of the first element in matr
matr[0] = 321;  // Change the first element
matr[1] = 322;  // Change the second element
printf("\nAddress of matrix element matr[1]: %p", &matr[0]);
printf("\nValue of the  matrix element  matr[1]; %d", matr[0]);
printf("\nAddress of matrix element matr[2]: %p", &matr[1]);
printf("\nValue of the matrix element  matr[2]: %d\n", matr[1]);
printf("\nValue of the pointer p: %p", p);
printf("\nThe value p points to: %d", *p);
printf("\nThe value that (p+1) points to  %d\n", *(p+1));
printf("\nAddress of pointer p : %p\n", &p);

Address of the matrix element matr[1]: 0xbfffef70
Value of the  matrix element  matr[1]; 321
Address of the matrix element matr[2]: 0xbfffef74
Value of the matrix element  matr[2]: 322
Value of the pointer: 0xbfffef70
The value pointer points at: 321
The value that (pointer+1) points at:  322
Address of the pointer variable : 0xbfffef6c

using namespace std;
#include <iostream>
int main(int argc, char *argv[])
{
  FILE *in_file, *out_file;
  if( argc < 3)  {
    printf("The programs has the following structure :\n");
    printf("write in the name of the input and output files \n");
    exit(0);
  }
  in_file = fopen( argv[1], "r");// returns pointer to the  input file
  if( in_file == NULL )  { // NULL means that the file is missing
    printf("Can't find the input file %s\n", argv[1]);
    exit(0);

 out_file = fopen( argv[2], "w"); // returns a pointer to the output file
 if( out_file == NULL )  {       // can't find the file
    printf("Can't find the output file%s\n", argv[2]);
    exit(0);
  }
  fclose(in_file);
  fclose(out_file);
  return 0;

#include <fstream>

// input and output file as global variable
ofstream ofile;
ifstream ifile;

int main(int argc, char* argv[])
{
  char *outfilename;
  //Read in output file, abort if there are too
  //few command-line arguments
  if( argc <= 1 ){
    cout << "Bad Usage: " << argv[0] <<
      " read also output file on same line" << endl;
    exit(1);
  }
  else{
    outfilename=argv[1];
  }
  ofile.open(outfilename);
  .....
  ofile.close();  // close output file

void output(double r_min , double r_max, int max_step,
            double *d)
{
int i;
ofile << "RESULTS:" << endl;
ofile << setiosflags(ios::showpoint | ios::uppercase);
ofile <<"R_min = " << setw(15) << setprecision(8) <<r_min <<endl;
ofile <<"R_max = " << setw(15) << setprecision(8) <<r_max <<endl;
ofile <<"Number of steps = " << setw(15) << max_step << endl;
ofile << "Five lowest eigenvalues:" << endl;
for(i = 0; i < 5; i++) {
    ofile << setw(15) << setprecision(8) << d[i] << endl;

}  // end of function output

int main(int argc, char* argv[])
{
  char *infilename;
  // Read in input file, abort if there are too
  // few command-line arguments
  if( argc <= 1 ){
    cout << "Bad Usage: " << argv[0] <<
      " read also input file on same line" << endl;
    exit(1);
  }
  else{
    infilename=argv[1];
  }
  ifile.open(infilename);
  ....
  ifile.close();  // close input file

const char* filename1 = "myfile";
ifstream ifile(filename1);
string filename2 = filename1 + ".out"
ofstream ofile(filename2);  // new output file
ofstream ofile(filename2, ios_base::app);  // append

//      Read something from the file:

double a; int b; char c[200];
ifile >> a >> b >> c;  // skips white space in between

//      Can test on success of reading:

if (!(ifile >> a >> b >> c)) ok = 0;

int main(int argc, char argv[]) {
int  a:
int *b;
a = 10;
b = new int[10];
for (i = 0; i < 10; i++) {
  b[i] = i;
}
func(a, b);
delete [] b;
return 0;
}

   int n; n =8;
   func(&n); /* &n is a pointer to n */
   ....
   void func(int *i)
   {
     *i = 10; /* n is changed to 10 */
     ....
   }

   int n; n =8;
   func(n); // just transfer n itself
   ....
   void func(int& i)
   {
     i = 10; // n is changed to 10
     ....
   }

SUBROUTINE coulomb_integral(np,lp,n,l,coulomb)
  USE effective_interaction_declar
  USE energy_variables
  USE wave_functions
  IMPLICIT NONE
  INTEGER, INTENT(IN)  :: n, l, np, lp
  INTEGER :: i
  REAL(KIND=8), INTENT(INOUT) :: coulomb
  REAL(KIND=8) :: z_rel, oscl_r, sum_coulomb
  ...

Matrix properties reminder.

$$ {\bf A} = \left( \begin{array}{cccc} a_{11} & a_{12} & a_{13} & a_{14} \\ a_{21} & a_{22} & a_{23} & a_{24} \\ a_{31} & a_{32} & a_{33} & a_{34} \\ a_{41} & a_{42} & a_{43} & a_{44} \end{array} \right)\qquad {\bf I} = \left( \begin{array}{cccc} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{array} \right) $$

The inverse of a matrix is defined by $$ {\bf A}^{-1} \cdot {\bf A} = I $$

Matrix Properties Reminder.

Relations	Name	matrix elements
\( A = A^{T} \)	symmetric	\( a_{ij} = a_{ji} \)
\( A = \left (A^{T} \right )^{-1} \)	real orthogonal	\( \sum_k a_{ik} a_{jk} = \sum_k a_{ki} a_{kj} = \delta_{ij} \)
\( A = A^{*} \)	real matrix	\( a_{ij} = a_{ij}^{*} \)
\( A = A^{\dagger} \)	hermitian	\( a_{ij} = a_{ji}^{*} \)
\( A = \left (A^{\dagger} \right )^{-1} \)	unitary	\( \sum_k a_{ik} a_{jk}^{*} = \sum_k a_{ki}^{*} a_{kj} = \delta_{ij} \)

Some Equivalent Statements.

For an \( N\times N \) matrix \( {\bf A} \) the following properties are all equivalent

• If the inverse of \( {\bf A} \) exists, \( {\bf A} \) is nonsingular.
• The equation \( {\bf Ax}=0 \) implies \( {\bf x}=0 \).
• The rows of \( {\bf A} \) form a basis of \( R^N \).
• The columns of \( {\bf A} \) form a basis of \( R^N \).
• \( {\bf A} \) is a product of elementary matrices.
• \( 0 \) is not eigenvalue of \( {\bf A} \).

Static.

We have an \( N\times N \) matrix A with \( N=100 \) In C/C++ this would be defined as

   int N = 100;
   double A[100][100];
   //   initialize all elements to zero
   for(i=0 ; i < N ; i++) {
      for(j=0 ; j < N ; j++) {
         A[i][j] = 0.0;

Note the way the matrix is organized, row-major order.

Row Major Order, Addition.

We have \( N\times N \) matrices A, B and C and we wish to evaluate \( A=B+C \). $$ {\bf A}= {\bf B}\pm{\bf C} \Longrightarrow a_{ij} = b_{ij}\pm c_{ij}, $$ In C/C++ this would be coded like

   for(i=0 ; i < N ; i++) {
      for(j=0 ; j < N ; j++) {
         a[i][j] = b[i][j]+c[i][j]

Row Major Order, Multiplication.

We have \( N\times N \) matrices A, B and C and we wish to evaluate \( A=BC \). $$ {\bf A}={\bf BC} \Longrightarrow a_{ij} = \sum_{k=1}^{n} b_{ik}c_{kj}, $$ In C/C++ this would be coded like

   for(i=0 ; i < N ; i++) {
      for(j=0 ; j < N ; j++) {
         for(k=0 ; k < N ; k++) {
            a[i][j]+=b[i][k]*c[k][j];

Column Major Order.

   ALLOCATE (a(N,N), b(N,N), c(N,N))
   DO j=1,  N
      DO i=1, N
         a(i,j)=b(i,j)+c(i,j)
      ENDDO
   ENDDO
   ...
   DEALLOCATE(a,b,c)

Fortran 90 writes the above statements in a much simpler way

   a=b+c

Multiplication

   a=MATMUL(b,c)

Fortran contains also the intrinsic functions TRANSPOSE and CONJUGATE.

Do it yourself.

int N;
double **  A;
A = new double*[N]
for ( i = 0; i < N; i++)
    A[i] = new double[N];

Always free space when you don't need an array anymore.

for ( i = 0; i < N; i++)
    delete[] A[i];
delete[] A;

#include <iostream>
#include <armadillo>

using namespace std;
using namespace arma;

int main(int argc, char** argv)
  {
  mat A = randu<mat>(5,5);
  mat B = randu<mat>(5,5);

  cout << A*B << endl;

  return 0;

c++ -O2 -o program.x program.cpp  -larmadillo -llapack -lblas

#include <iostream>
#include "armadillo"
using namespace arma;
using namespace std;

int main(int argc, char** argv)
  {
  // directly specify the matrix size (elements are uninitialised)
  mat A(2,3);
  // .n_rows = number of rows    (read only)
  // .n_cols = number of columns (read only)
  cout << "A.n_rows = " << A.n_rows << endl;
  cout << "A.n_cols = " << A.n_cols << endl;
  // directly access an element (indexing starts at 0)
  A(1,2) = 456.0;
  A.print("A:");
  // scalars are treated as a 1x1 matrix,
  // hence the code below will set A to have a size of 1x1
  A = 5.0;
  A.print("A:");
  // if you want a matrix with all elements set to a particular value
  // the .fill() member function can be used
  A.set_size(3,3);
  A.fill(5.0);  A.print("A:");

  mat B;

  // endr indicates "end of row"
  B << 0.555950 << 0.274690 << 0.540605 << 0.798938 << endr
    << 0.108929 << 0.830123 << 0.891726 << 0.895283 << endr
    << 0.948014 << 0.973234 << 0.216504 << 0.883152 << endr
    << 0.023787 << 0.675382 << 0.231751 << 0.450332 << endr;

  // print to the cout stream
  // with an optional string before the contents of the matrix
  B.print("B:");

  // the << operator can also be used to print the matrix
  // to an arbitrary stream (cout in this case)
  cout << "B:" << endl << B << endl;
  // save to disk
  B.save("B.txt", raw_ascii);
  // load from disk
  mat C;
  C.load("B.txt");
  C += 2.0 * B;
  C.print("C:");

  // submatrix types:
  //
  // .submat(first_row, first_column, last_row, last_column)
  // .row(row_number)
  // .col(column_number)
  // .cols(first_column, last_column)
  // .rows(first_row, last_row)

  cout << "C.submat(0,0,3,1) =" << endl;
  cout << C.submat(0,0,3,1) << endl;

  // generate the identity matrix
  mat D = eye<mat>(4,4);

  D.submat(0,0,3,1) = C.cols(1,2);
  D.print("D:");

  // transpose
  cout << "trans(B) =" << endl;
  cout << trans(B) << endl;

  // maximum from each column (traverse along rows)
  cout << "max(B) =" << endl;
  cout << max(B) << endl;

  // maximum from each row (traverse along columns)
  cout << "max(B,1) =" << endl;
  cout << max(B,1) << endl;
  // maximum value in B
  cout << "max(max(B)) = " << max(max(B)) << endl;
  // sum of each column (traverse along rows)
  cout << "sum(B) =" << endl;
  cout << sum(B) << endl;
  // sum of each row (traverse along columns)
  cout << "sum(B,1) =" << endl;
  cout << sum(B,1) << endl;
  // sum of all elements
  cout << "sum(sum(B)) = " << sum(sum(B)) << endl;
  cout << "accu(B)     = " << accu(B) << endl;
  // trace = sum along diagonal
  cout << "trace(B)    = " << trace(B) << endl;
  // random matrix -- values are uniformly distributed in the [0,1] interval
  mat E = randu<mat>(4,4);
  E.print("E:");

  // row vectors are treated like a matrix with one row
  rowvec r;
  r << 0.59499 << 0.88807 << 0.88532 << 0.19968;
  r.print("r:");

  // column vectors are treated like a matrix with one column
  colvec q;
  q << 0.81114 << 0.06256 << 0.95989 << 0.73628;
  q.print("q:");

  // dot or inner product
  cout << "as_scalar(r*q) = " << as_scalar(r*q) << endl;

    // outer product
  cout << "q*r =" << endl;
  cout << q*r << endl;


  // sum of three matrices (no temporary matrices are created)
  mat F = B + C + D;
  F.print("F:");

    return 0;

#include <iostream>
#include "armadillo"
using namespace arma;
using namespace std;

int main(int argc, char** argv)
  {
  cout << "Armadillo version: " << arma_version::as_string() << endl;

  mat A;

  A << 0.165300 << 0.454037 << 0.995795 << 0.124098 << 0.047084 << endr
    << 0.688782 << 0.036549 << 0.552848 << 0.937664 << 0.866401 << endr
    << 0.348740 << 0.479388 << 0.506228 << 0.145673 << 0.491547 << endr
    << 0.148678 << 0.682258 << 0.571154 << 0.874724 << 0.444632 << endr
    << 0.245726 << 0.595218 << 0.409327 << 0.367827 << 0.385736 << endr;

  A.print("A =");

  // determinant
  cout << "det(A) = " << det(A) << endl;

  // inverse
  cout << "inv(A) = " << endl << inv(A) << endl;
  double k = 1.23;

  mat    B = randu<mat>(5,5);
  mat    C = randu<mat>(5,5);

  rowvec r = randu<rowvec>(5);
  colvec q = randu<colvec>(5);


  // examples of some expressions
  // for which optimised implementations exist
  // optimised implementation of a trinary expression
  // that results in a scalar
  cout << "as_scalar( r*inv(diagmat(B))*q ) = ";
  cout << as_scalar( r*inv(diagmat(B))*q ) << endl;

  // example of an expression which is optimised
  // as a call to the dgemm() function in BLAS:
  cout << "k*trans(B)*C = " << endl << k*trans(B)*C;

    return 0;

   btemp = b[1];
   u[1] = f[1]/btemp;
   for(i=2 ; i <= n ; i++) {
      temp[i] = c[i-1]/btemp;
      btemp = b[i]-a[i]*temp[i];
      u[i] = (f[i] - a[i]*u[i-1])/btemp;

   for(i=n-1 ; i >= 1 ; i--) {
      u[i] -= temp[i+1]*u[i+1];

//  Simple matrix inversion example
#include <iostream>
#include <new>
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <cstring>
#include "lib.h"

using namespace std;

/* function declarations */

void inverse(double **, int);

void inverse(double **a, int n)
{
  int          i,j, *indx;
  double       d, *col, **y;
  // allocate space in memory
  indx = new int[n];
  col  = new double[n];
  y    = (double **) matrix(n, n, sizeof(double));
  ludcmp(a, n, indx, &d);   // LU decompose  a[][]
  printf("\n\nLU form of matrix of a[][]:\n");
  for(i = 0; i < n; i++) {
    printf("\n");
    for(j = 0; j < n; j++) {
      printf(" a[%2d][%2d] = %12.4E",i, j, a[i][j]);

  // find inverse of a[][] by columns
  for(j = 0; j < n; j++) {
    // initialize right-side of linear equations
    for(i = 0; i < n; i++) col[i] = 0.0;
    col[j] = 1.0;
    lubksb(a, n, indx, col);
    // save result in y[][]
    for(i = 0; i < n; i++) y[i][j] = col[i];
  }   //j-loop over columns
  // return the inverse matrix in a[][]
  for(i = 0; i < n; i++) {
    for(j = 0; j < n; j++) a[i][j] = y[i][j];

  free_matrix((void **) y);     // release local memory
  delete [] col;
  delete []indx;
}  // End: function inverse()

PROGRAM matrix
  USE constants
  USE F90library
  IMPLICIT NONE
  !      The definition of the matrix, using dynamic allocation
  REAL(DP), ALLOCATABLE, DIMENSION(:,:) :: a, ainv, unity
  !      the determinant
  REAL(DP) :: d
  !      The size of the matrix
  INTEGER :: n
  ....
  !      Allocate now place in heap for a
  ALLOCATE ( a(n,n), ainv(n,n), unity(n,n) )

  WRITE(6,*) ' The matrix before inversion'
  WRITE(6,'(3F12.6)') a
  ainv=a
  CALL matinv (ainv, n, d)
  ....
  !      get the unity matrix
  unity=MATMUL(ainv,a)
  WRITE(6,*) ' The unity matrix'
  WRITE(6,'(3F12.6)') unity
  !      deallocate all arrays
  DEALLOCATE (a, ainv, unity)
END PROGRAM matrix

Linear Algebra and project 1.

• Discussion of Project 1
• Object orientation and unit testing (see detailed instruction on webpage)

   for(i=0 ; i < n ; i++) {
      for(j=0 ; j < n ; j++) {
         a[i][j]=b[i][j]+c[i][j]

// Program to calculate addition and multiplication of two complex numbers
using namespace std;
#include <iostream>
#include <cmath>
#include <complex>
int main()
{
  complex<double> x(6.1,8.2), y(0.5,1.3);
  // write out x+y
  cout << x + y << x*y  << endl;
  return 0;

#ifndef Complex_H
#define Complex_H
//   various include statements and definitions
#include <iostream>          // Standard ANSI-C++ include files
#include <new>
#include ....

class Complex
{...
definition of variables and their character
};
//   declarations of various functions used by the class
...
#endif

#include "Complex.h"
...  other include and declarations
int main ()
{
  Complex a(0.1,1.3);    // we declare a complex variable a
  Complex b(3.0), c(5.0,-2.3);  // we declare  complex variables b and c
  Complex d = b;         //  we declare  a new complex variable d
  cout << "d=" << d << ", a=" << a << ", b=" << b << endl;
  d = a*c + b/a;  //   we add, multiply and divide two complex numbers
  cout << "Re(d)=" << d.Re() << ", Im(d)=" << d.Im() << endl;  // write out of the real and imaginary parts

class Complex
{
private:
   double re, im; // real and imaginary part
public:
   Complex ();                              // Complex c;
   Complex (double re, double im = 0.0); // Definition of a complex variable;
   Complex (const Complex& c);              // Usage: Complex c(a);   // equate two complex variables
   Complex& operator= (const Complex& c); // c = a;   //  equate two complex variables, same as previous
....

  ~Complex () {}                        // destructor
   double   Re () const;        // double real_part = a.Re();
   double   Im () const;        // double imag_part = a.Im();
   double   abs () const;       // double m = a.abs(); // modulus
   friend Complex operator+ (const Complex&  a, const Complex& b);
   friend Complex operator- (const Complex&  a, const Complex& b);
   friend Complex operator* (const Complex&  a, const Complex& b);
   friend Complex operator/ (const Complex&  a, const Complex& b);
};

Complex:: Complex () { re = im = 0.0; }

Complex:: Complex () {}

Complex:: Complex (double re_a, double im_a) {
    re = re_a; im = im_a; }

double Complex:: Re () const { return re; }} //  getting the real part
double Complex:: Im () const { return im; }  //   and the imaginary part

// const arguments (in functions) cannot be changed:
void myfunc (const Complex& c)
{ c.re = 0.2; /* ILLEGAL!! compiler error... */  }

double Complex:: abs ()  { return sqrt(re*re + im*im);}

double myabs (const Complex& c)
{ return c.abs(); }   // Not ok because c.abs() is not a const func.

double Complex:: abs () const { return sqrt(re*re + im*im); }

Complex& Complex:: operator= (const Complex& c)
{
   re = c.re;
   im = c.im;
   return *this;
}

Complex:: Complex (const Complex& c)
{ *this = c; }

inline Complex operator+ (const Complex& a, const Complex& b)
{ return Complex (a.re + b.re, a.im + b.im); }

c.re = a.re + b.re;
c.im = a.im + b.im;

class Complex
{
   ...
   friend Complex operator+ (const Complex& a, const Complex& b);
   ...
};

inline Complex operator+ (const Complex& a, const Complex& b)
{ return Complex (a.re + b.re, a.im + b.im); }

inline Complex operator+ (const Complex& a, const Complex& b)
{ return Complex (a.Re() + b.Re(), a.Im() + b.Im()); }

inline Complex operator* (const Complex& a, const Complex& b)
{
  return Complex(a.re*b.re - a.im*b.im, a.im*b.re + a.re*b.im);

inline Complex:: Complex () { re = im = 0.0; }
inline Complex:: Complex (double re_, double im_)
{ ... }
inline Complex:: Complex (const Complex& c)
{ ... }
inline Complex:: operator= (const Complex& c)
{ ... }

// e, c, d are complex
e = c*d;
// first compiler translation:
e.operator= (operator* (c,d));
// result of nested inline functions
// operator=, operator*, Complex(double,double=0):
e.re = c.re*d.re - c.im*d.im;
e.im = c.im*d.re + c.re*d.im;

ostream& operator<< (ostream& o, const Complex& c)
{ o << "(" << c.Re() << "," << c.Im() << ") "; return o;}

template<class T>
class Complex
{
private:
   T re, im; // real and imaginary part
public:
   Complex ();                              // Complex c;
   Complex (T re, T im = 0); // Definition of a complex variable;
   Complex (const Complex& c);              // Usage: Complex c(a);   // equate two complex variables
   Complex& operator= (const Complex& c); // c = a;   //  equate two complex variables, same as previous

  ~Complex () {}                        // destructor
   T   Re () const;        // T real_part = a.Re();
   T   Im () const;        // T imag_part = a.Im();
   T   abs () const;       // T m = a.abs(); // modulus
   friend Complex operator+ (const Complex&  a, const Complex& b);
   friend Complex operator- (const Complex&  a, const Complex& b);
   friend Complex operator* (const Complex&  a, const Complex& b);
   friend Complex operator/ (const Complex&  a, const Complex& b);
};

Complex<double> a(10.0,5.1);
Complex<int> b(1,0);

template<class T>
Complex<T>:: Complex (T re_a, T im_a)
{ re = re_a; im = im_a; }

// Not all declarations here
// Class to compute the square of a number
template<class T>
class Squared{
  public:
    // Default constructor, not used here
    Squared(){}

    // Overload the function operator()
    T operator()(T x){return x*x;}

};

#include <iostream>
#include "squared.h"
using namespace std;

int main(){
  Squared<double> s;
  cout << s(3) << endl;

Eigenvalue problems and project 2.

• Discussion of Jacobi's algorithm, chapter 7.1-7.2
• Presentation of project 2.

Reading assignment this week: chapters 7.1-7.4 of lecture notes

Let us consider the matrix \( {\bf A} \) of dimension \( n \). The eigenvalues of \( {\bf A} \) are defined through the matrix equation $$ {\bf A}{\bf x}^{(\nu)} = \lambda^{(\nu)}{\bf x}^{(\nu)}, $$ where \( \lambda^{(\nu)} \) are the eigenvalues and \( {\bf x}^{(\nu)} \) the corresponding eigenvectors. Unless otherwise stated, when we use the wording eigenvector we mean the right eigenvector. The left eigenvalue problem is defined as $$ {\bf x}^{(\nu)}_L{\bf A} = \lambda^{(\nu)}{\bf x}^{(\nu)}_L $$ The above right eigenvector problem is equivalent to a set of \( n \) equations with \( n \) unknowns \( x_i \).

The eigenvalue problem can be rewritten as $$ \left( {\bf A}-\lambda^{(\nu)} {\bf I} \right) {\bf x}^{(\nu)} = 0, $$ with \( {\bf I} \) being the unity matrix. This equation provides a solution to the problem if and only if the determinant is zero, namely $$ \left| {\bf A}-\lambda^{(\nu)}{\bf I}\right| = 0, $$ which in turn means that the determinant is a polynomial of degree \( n \) in \( \lambda \) and in general we will have \( n \) distinct zeros.

The eigenvalues of a matrix \( {\bf A}\in {\mathbb{C}}^{n\times n} \) are thus the \( n \) roots of its characteristic polynomial $$ P(\lambda) = det(\lambda{\bf I}-{\bf A}), $$ or $$ P(\lambda)= \prod_{i=1}^{n}\left(\lambda_i-\lambda\right). $$ The set of these roots is called the spectrum and is denoted as \( \lambda({\bf A}) \). If \( \lambda({\bf A})=\left\{\lambda_1,\lambda_2,\dots ,\lambda_n\right\} \) then we have $$ det({\bf A})= \lambda_1\lambda_2\dots\lambda_n, $$ and if we define the trace of \( {\bf A} \) as $$ Tr({\bf A})=\sum_{i=1}^n a_{ii}$$ then $$ Tr({\bf A})=\lambda_1+\lambda_2+\dots+\lambda_n. $$

The Abel-Ruffini theorem (also known as Abel's impossibility theorem) states that there is no general solution in radicals to polynomial equations of degree five or higher.

The content of this theorem is frequently misunderstood. It does not assert that higher-degree polynomial equations are unsolvable. In fact, if the polynomial has real or complex coefficients, and we allow complex solutions, then every polynomial equation has solutions; this is the fundamental theorem of algebra. Although these solutions cannot always be computed exactly with radicals, they can be computed to any desired degree of accuracy using numerical methods such as the Newton-Raphson method or Laguerre method, and in this way they are no different from solutions to polynomial equations of the second, third, or fourth degrees.

The theorem only concerns the form that such a solution must take. The content of the theorem is that the solution of a higher-degree equation cannot in all cases be expressed in terms of the polynomial coefficients with a finite number of operations of addition, subtraction, multiplication, division and root extraction. Some polynomials of arbitrary degree, of which the simplest nontrivial example is the monomial equation \( ax^n = b \), are always solvable with a radical.

The Abel-Ruffini theorem says that there are some fifth-degree equations whose solution cannot be so expressed. The equation \( x^5 - x + 1 = 0 \) is an example. Some other fifth degree equations can be solved by radicals, for example \( x^5 - x^4 - x + 1 = 0 \). The precise criterion that distinguishes between those equations that can be solved by radicals and those that cannot was given by Galois and is now part of Galois theory: a polynomial equation can be solved by radicals if and only if its Galois group is a solvable group.

Today, in the modern algebraic context, we say that second, third and fourth degree polynomial equations can always be solved by radicals because the symmetric groups \( S_2, S_3 \) and \( S_4 \) are solvable groups, whereas \( S_n \) is not solvable for \( n \ge 5 \).

In the present discussion we assume that our matrix is real and symmetric, that is \( {\bf A}\in {\mathbb{R}}^{n\times n} \). The matrix \( {\bf A} \) has \( n \) eigenvalues \( \lambda_1\dots \lambda_n \) (distinct or not). Let \( {\bf D} \) be the diagonal matrix with the eigenvalues on the diagonal $$ {\bf D}= \left( \begin{array}{ccccccc} \lambda_1 & 0 & 0 & 0 & \dots &0 & 0 \\ 0 & \lambda_2 & 0 & 0 & \dots &0 &0 \\ 0 & 0 & \lambda_3 & 0 &0 &\dots & 0\\ \dots & \dots & \dots & \dots &\dots &\dots & \dots\\ 0 & \dots & \dots & \dots &\dots &\lambda_{n-1} & \\ 0 & \dots & \dots & \dots &\dots &0 & \lambda_n \end{array} \right). $$ If \( {\bf A} \) is real and symmetric then there exists a real orthogonal matrix \( {\bf S} \) such that $$ {\bf S}^T {\bf A}{\bf S}= \mathrm{diag}(\lambda_1,\lambda_2,\dots ,\lambda_n), $$ and for \( j=1:n \) we have \( {\bf A}{\bf S}(:,j) = \lambda_j {\bf S}(:,j) \).

To obtain the eigenvalues of \( {\bf A}\in {\mathbb{R}}^{n\times n} \), the strategy is to perform a series of similarity transformations on the original matrix \( {\bf A} \), in order to reduce it either into a diagonal form as above or into a tridiagonal form.

We say that a matrix \( {\bf B} \) is a similarity transform of \( {\bf A} \) if $$ {\bf B}= {\bf S}^T {\bf A}{\bf S}, \hspace{1cm} \mathrm{where} \hspace{1cm} {\bf S}^T{\bf S}={\bf S}^{-1}{\bf S} ={\bf I}. $$ The importance of a similarity transformation lies in the fact that the resulting matrix has the same eigenvalues, but the eigenvectors are in general different.

To prove this we start with the eigenvalue problem and a similarity transformed matrix \( {\bf B} \). $$ {\bf A}{\bf x}=\lambda{\bf x} \hspace{1cm} \mathrm{and}\hspace{1cm} {\bf B}= {\bf S}^T {\bf A}{\bf S}. $$ We multiply the first equation on the left by \( {\bf S}^T \) and insert \( {\bf S}^{T}{\bf S} = {\bf I} \) between \( {\bf A} \) and \( {\bf x} \). Then we get $$ \begin{equation} ({\bf S}^T{\bf A}{\bf S})({\bf S}^T{\bf x})=\lambda{\bf S}^T{\bf x} , \end{equation} $$ which is the same as $$ {\bf B} \left ( {\bf S}^T{\bf x} \right ) = \lambda \left ({\bf S}^T{\bf x}\right ). $$ The variable \( \lambda \) is an eigenvalue of \( {\bf B} \) as well, but with eigenvector \( {\bf S}^T{\bf x} \).

The basic philosophy is to

• Either apply subsequent similarity transformations (direct method) so that

$$ \begin{equation} {\bf S}_N^T\dots {\bf S}_1^T{\bf A}{\bf S}_1\dots {\bf S}_N={\bf D} , \end{equation} $$

• Or apply subsequent similarity transformations so that \( {\bf A} \) becomes tridiagonal (Householder) or upper/lower triangular (the QR method to be discussed later).
• Thereafter, techniques for obtaining eigenvalues from tridiagonal matrices can be used.
• Or use so-called power methods
• Or use iterative methods (Krylov, Lanczos, Arnoldi). These methods are popular for huge matrix problems.

In project 1 we rewrote our original differential equation in terms of a discretized equation with approximations to the derivatives as $$ -\frac{u_{i+1} -2u_i +u_{i-i}}{h^2}=f(x_i,u(x_i)), $$ with \( i=1,2,\dots, n \). We need to add to this system the two boundary conditions \( u(a) =u_0 \) and \( u(b) = u_{n+1} \). If we define a matrix $$ {\bf A} = \frac{1}{h^2}\left(\begin{array}{cccccc} 2 & -1 & & & & \\ -1 & 2 & -1 & & & \\ & -1 & 2 & -1 & & \\ & \dots & \dots &\dots &\dots & \dots \\ & & &-1 &2& -1 \\ & & & &-1 & 2 \\ \end{array} \right) $$ and the corresponding vectors \( {\bf u} = (u_1, u_2, \dots,u_n)^T \) and \( {\bf f}({\bf u}) = f(x_1,x_2,\dots, x_n,u_1, u_2, \dots,u_n)^T \) we can rewrite the differential equation including the boundary conditions as a system of linear equations with a large number of unknowns $$ {\bf A}{\bf u} = {\bf f}({\bf u}). $$

We are first interested in the solution of the radial part of Schroedinger's equation for one electron. This equation reads $$ -\frac{\hbar^2}{2 m} \left ( \frac{1}{r^2} \frac{d}{dr} r^2 \frac{d}{dr} - \frac{l (l + 1)}{r^2} \right )R(r) + V(r) R(r) = E R(r). $$ In our case \( V(r) \) is the harmonic oscillator potential \( (1/2)kr^2 \) with \( k=m\omega^2 \) and \( E \) is the energy of the harmonic oscillator in three dimensions. The oscillator frequency is \( \omega \) and the energies are $$ E_{nl}= \hbar \omega \left(2n+l+\frac{3}{2}\right), $$ with \( n=0,1,2,\dots \) and \( l=0,1,2,\dots \).

Since we have made a transformation to spherical coordinates it means that \( r\in [0,\infty) \). The quantum number \( l \) is the orbital momentum of the electron. Then we substitute \( R(r) = (1/r) u(r) \) and obtain $$ -\frac{\hbar^2}{2 m} \frac{d^2}{dr^2} u(r) + \left ( V(r) + \frac{l (l + 1)}{r^2}\frac{\hbar^2}{2 m} \right ) u(r) = E u(r) . $$ The boundary conditions are \( u(0)=0 \) and \( u(\infty)=0 \).

We introduce a dimensionless variable \( \rho = (1/\alpha) r \) where \( \alpha \) is a constant with dimension length and get $$ -\frac{\hbar^2}{2 m \alpha^2} \frac{d^2}{d\rho^2} u(\rho) + \left ( V(\rho) + \frac{l (l + 1)}{\rho^2} \frac{\hbar^2}{2 m\alpha^2} \right ) u(\rho) = E u(\rho) . $$ In project 2 we choose \( l=0 \). Inserting \( V(\rho) = (1/2) k \alpha^2\rho^2 \) we end up with $$ -\frac{\hbar^2}{2 m \alpha^2} \frac{d^2}{d\rho^2} u(\rho) + \frac{k}{2} \alpha^2\rho^2u(\rho) = E u(\rho) . $$ We multiply thereafter with \( 2m\alpha^2/\hbar^2 \) on both sides and obtain $$ -\frac{d^2}{d\rho^2} u(\rho) + \frac{mk}{\hbar^2} \alpha^4\rho^2u(\rho) = \frac{2m\alpha^2}{\hbar^2}E u(\rho) . $$

We have thus $$ -\frac{d^2}{d\rho^2} u(\rho) + \frac{mk}{\hbar^2} \alpha^4\rho^2u(\rho) = \frac{2m\alpha^2}{\hbar^2}E u(\rho) . $$ The constant \( \alpha \) can now be fixed so that $$ \frac{mk}{\hbar^2} \alpha^4 = 1, $$ or $$ \alpha = \left(\frac{\hbar^2}{mk}\right)^{1/4}. $$ Defining $$ \lambda = \frac{2m\alpha^2}{\hbar^2}E, $$ we can rewrite Schr\"odinger's equation as $$ -\frac{d^2}{d\rho^2} u(\rho) + \rho^2u(\rho) = \lambda u(\rho) . $$ This is the first equation to solve numerically. In three dimensions the eigenvalues for \( l=0 \) are \( \lambda_0=3,\lambda_1=7,\lambda_2=11,\dots . \)

We use the by now standard expression for the second derivative of a function \( u \) $$ \begin{equation} u''=\frac{u(\rho+h) -2u(\rho) +u(\rho-h)}{h^2} +O(h^2), \label{eq:diffoperation} \end{equation} $$ where \( h \) is our step. Next we define minimum and maximum values for the variable \( \rho \), \( \rho_{\mathrm{min}}=0 \) and \( \rho_{\mathrm{max}} \), respectively. You need to check your results for the energies against different values \( \rho_{\mathrm{max}} \), since we cannot set \( \rho_{\mathrm{max}}=\infty \).

With a given number of steps, \( n_{\mathrm{step}} \), we then define the step \( h \) as $$ h=\frac{\rho_{\mathrm{max}}-\rho_{\mathrm{min}} }{n_{\mathrm{step}}}. $$ Define an arbitrary value of \( \rho \) as $$ \rho_i= \rho_{\mathrm{min}} + ih \hspace{1cm} i=0,1,2,\dots , n_{\mathrm{step}} $$ we can rewrite the Schr\"odinger equation for \( \rho_i \) as $$ -\frac{u(\rho_i+h) -2u(\rho_i) +u(\rho_i-h)}{h^2}+\rho_i^2u(\rho_i) = \lambda u(\rho_i), $$ or in a more compact way $$ -\frac{u_{i+1} -2u_i +u_{i-1}}{h^2}+\rho_i^2u_i=-\frac{u_{i+1} -2u_i +u_{i-1} }{h^2}+V_iu_i = \lambda u_i, $$ where \( V_i=\rho_i^2 \) is the harmonic oscillator potential.

Define first the diagonal matrix element $$ d_i=\frac{2}{h^2}+V_i, $$ and the non-diagonal matrix element $$ e_i=-\frac{1}{h^2}. $$ In this case the non-diagonal matrix elements are given by a mere constant. All non-diagonal matrix elements are equal.

With these definitions the Schroedinger equation takes the following form $$ d_iu_i+e_{i-1}u_{i-1}+e_{i+1}u_{i+1} = \lambda u_i, $$ where \( u_i \) is unknown. We can write the latter equation as a matrix eigenvalue problem $$ \begin{equation} \left( \begin{array}{ccccccc} d_1 & e_1 & 0 & 0 & \dots &0 & 0 \\ e_1 & d_2 & e_2 & 0 & \dots &0 &0 \\ 0 & e_2 & d_3 & e_3 &0 &\dots & 0\\ \dots & \dots & \dots & \dots &\dots &\dots & \dots\\ 0 & \dots & \dots & \dots &\dots &d_{n_{\mathrm{step}}-2} & e_{n_{\mathrm{step}}-1}\\ 0 & \dots & \dots & \dots &\dots &e_{n_{\mathrm{step}}-1} & d_{n_{\mathrm{step}}-1} \end{array} \right) \left( \begin{array}{c} u_{1} \\ u_{2} \\ \dots\\ \dots\\ \dots\\ u_{n_{\mathrm{step}}-1} \end{array} \right)=\lambda \left( \begin{array}{c} u_{1} \\ u_{2} \\ \dots\\ \dots\\ \dots\\ u_{n_{\mathrm{step}}-1} \end{array} \right) \label{eq:sematrix} \end{equation} $$ or if we wish to be more detailed, we can write the tridiagonal matrix as $$ \begin{equation} \left( \begin{array}{ccccccc} \frac{2}{h^2}+V_1 & -\frac{1}{h^2} & 0 & 0 & \dots &0 & 0 \\ -\frac{1}{h^2} & \frac{2}{h^2}+V_2 & -\frac{1}{h^2} & 0 & \dots &0 &0 \\ 0 & -\frac{1}{h^2} & \frac{2}{h^2}+V_3 & -\frac{1}{h^2} &0 &\dots & 0\\ \dots & \dots & \dots & \dots &\dots &\dots & \dots\\ 0 & \dots & \dots & \dots &\dots &\frac{2}{h^2}+V_{n_{\mathrm{step}}-2} & -\frac{1}{h^2}\\ 0 & \dots & \dots & \dots &\dots &-\frac{1}{h^2} & \frac{2}{h^2}+V_{n_{\mathrm{step}}-1} \end{array} \right) \label{eq:matrixse} \end{equation} $$ Recall that the solutions are known via the boundary conditions at \( i=n_{\mathrm{step}} \) and at the other end point, that is for \( \rho_0 \). The solution is zero in both cases.

We are going to study two electrons in a harmonic oscillator well which also interact via a repulsive Coulomb interaction. Let us start with the single-electron equation written as $$ -\frac{\hbar^2}{2 m} \frac{d^2}{dr^2} u(r) + \frac{1}{2}k r^2u(r) = E^{(1)} u(r), $$ where \( E^{(1)} \) stands for the energy with one electron only. For two electrons with no repulsive Coulomb interaction, we have the following Schroedinger equation $$ \left( -\frac{\hbar^2}{2 m} \frac{d^2}{dr_1^2} -\frac{\hbar^2}{2 m} \frac{d^2}{dr_2^2}+ \frac{1}{2}k r_1^2+ \frac{1}{2}k r_2^2\right)u(r_1,r_2) = E^{(2)} u(r_1,r_2) . $$

Note that we deal with a two-electron wave function \( u(r_1,r_2) \) and two-electron energy \( E^{(2)} \).

With no interaction this can be written out as the product of two single-electron wave functions, that is we have a solution on closed form.

We introduce the relative coordinate \( {\bf r} = {\bf r}_1-{\bf r}_2 \) and the center-of-mass coordinate \( {\bf R} = 1/2({\bf r}_1+{\bf r}_2) \). With these new coordinates, the radial Schr\"odinger equation reads $$ \left( -\frac{\hbar^2}{m} \frac{d^2}{dr^2} -\frac{\hbar^2}{4 m} \frac{d^2}{dR^2}+ \frac{1}{4} k r^2+ kR^2\right)u(r,R) = E^{(2)} u(r,R). $$

The equations for \( r \) and \( R \) can be separated via the ansatz for the wave function \( u(r,R) = \psi(r)\phi(R) \) and the energy is given by the sum of the relative energy \( E_r \) and the center-of-mass energy \( E_R \), that is $$ E^{(2)}=E_r+E_R. $$ We add then the repulsive Coulomb interaction between two electrons, namely a term $$ V(r_1,r_2) = \frac{\beta e^2}{|{\bf r}_1-{\bf r}_2|}=\frac{\beta e^2}{r}, $$ with \( \beta e^2=1.44 \) eVnm.

Adding this term, the \( r \)-dependent Schroedinger equation becomes $$ \left( -\frac{\hbar^2}{m} \frac{d^2}{dr^2}+ \frac{1}{4}k r^2+\frac{\beta e^2}{r}\right)\psi(r) = E_r \psi(r). $$ This equation is similar to the one we had previously in parts (a) and (b) and we introduce again a dimensionless variable \( \rho = r/\alpha \). Repeating the same steps, we arrive at $$ -\frac{d^2}{d\rho^2} \psi(\rho) + \frac{mk}{4\hbar^2} \alpha^4\rho^2\psi(\rho)+\frac{m\alpha \beta e^2}{\rho\hbar^2}\psi(\rho) = \frac{m\alpha^2}{\hbar^2}E_r \psi(\rho) . $$

We want to manipulate this equation further to make it as similar to that in (a) as possible. We define a 'frequency' $$ \omega_r^2=\frac{1}{4}\frac{mk}{\hbar^2} \alpha^4, $$ and fix the constant \( \alpha \) by requiring $$ \frac{m\alpha \beta e^2}{\hbar^2}=1 $$ or $$ \alpha = \frac{\hbar^2}{m\beta e^2}. $$

Defining $$ \lambda = \frac{m\alpha^2}{\hbar^2}E, $$ we can rewrite Schroedinger's equation as $$ -\frac{d^2}{d\rho^2} \psi(\rho) + \omega_r^2\rho^2\psi(\rho) +\frac{1}{\rho}\psi(\rho) = \lambda \psi(\rho). $$

We treat \( \omega_r \) as a parameter which reflects the strength of the oscillator potential.

Here we will study the cases \( \omega_r = 0.01 \), \( \omega_r = 0.5 \), \( \omega_r =1 \), and \( \omega_r = 5 \) for the ground state only, that is the lowest-lying state.

With no repulsive Coulomb interaction you should get a result which corresponds to the relative energy of a non-interacting system. Make sure your results are stable as functions of \( \rho_{\mathrm{max}} \) and the number of steps.

We are only interested in the ground state with \( l=0 \). We omit the center-of-mass energy.

For specific oscillator frequencies, the above equation has analytic answers, see the article by M. Taut, Phys. Rev. A 48, 3561 - 3566 (1993). The article can be retrieved from the following web address http://prola.aps.org/abstract/PRA/v48/i5/p3561_1.

The general overview.

One speaks normally of two main approaches to solving the eigenvalue problem.

• The first is the formal method, involving determinants and the characteristic polynomial. This proves how many eigenvalues there are, and is the way most of you learned about how to solve the eigenvalue problem, but for matrices of dimensions greater than 2 or 3, it is rather impractical.
• The other general approach is to use similarity or unitary tranformations to reduce a matrix to diagonal form. This is normally done in two steps: first reduce to for example a tridiagonal form, and then to diagonal form. The main algorithms we will discuss in detail, Jacobi's and Householder's (so-called direct method) and Lanczos algorithms (an iterative method), follow this

methodology.

Direct or non-iterative methods require for matrices of dimensionality \( n\times n \) typically \( O(n^3) \) operations. These methods are normally called standard methods and are used for dimensionalities \( n \sim 10^5 \) or smaller. A brief historical overview

Year	\( n \)
1950	\( n=20 \)	(Wilkinson)
1965	\( n=200 \)	(Forsythe et al.)
1980	\( n=2000 \)	Linpack
1995	\( n=20000 \)	Lapack
2012	\( n\sim 10^5 \)	Lapack

shows that in the course of 60 years the dimension that direct diagonalization methods can handle has increased by almost a factor of \( 10^4 \). However, it pales beside the progress achieved by computer hardware, from flops to petaflops, a factor of almost \( 10^{15} \). We see clearly played out in history the \( O(n^3) \) bottleneck of direct matrix algorithms.

Sloppily speaking, when \( n\sim 10^4 \) is cubed we have \( O(10^{12}) \) operations, which is smaller than the \( 10^{15} \) increase in flops.

If the matrix to diagonalize is large and sparse, direct methods simply become impractical, also because many of the direct methods tend to destroy sparsity. As a result large dense matrices may arise during the diagonalization procedure. The idea behind iterative methods is to project the $n-$dimensional problem in smaller spaces, so-called Krylov subspaces. Given a matrix \( {\bf A} \) and a vector \( {\bf v} \), the associated Krylov sequences of vectors (and thereby subspaces) \( {\bf v} \), \( {\bf A}{\bf v} \), \( {\bf A}^2{\bf v} \), \( {\bf A}^3{\bf v},\dots \), represent successively larger Krylov subspaces.

Matrix	\( {\bf A}{\bf x}={\bf b} \)	\( {\bf A}{\bf x}=\lambda{\bf x} \)
\( {\bf A}={\bf A}^* \)	Conjugate gradient	Lanczos
\( {\bf A}\ne {\bf A}^* \)	GMRES etc	Arnoldi

The Numerical Recipes codes have been rewritten in Fortran 90/95 and C/C++ by us. The original source codes are taken from the widely used software package LAPACK, which follows two other popular packages developed in the 1970s, namely EISPACK and LINPACK.

• LINPACK: package for linear equations and least square problems.
• LAPACK:package for solving symmetric, unsymmetric and generalized eigenvalue problems. From LAPACK's website http://www.netlib.org it is possible to download for free all source codes from this library. Both C/C++ and Fortran versions are available.
• BLAS (I, II and III): (Basic Linear Algebra Subprograms) are routines that provide standard building blocks for performing basic vector and matrix operations. Blas I is vector operations, II vector-matrix operations and III matrix-matrix operations.

Consider an example of an (\( n\times n \)) orthogonal transformation matrix $$ {\bf S}= \left( \begin{array}{cccccccc} 1 & 0 & \dots & 0 & 0 & \dots & 0 & 0 \\ 0 & 1 & \dots & 0 & 0 & \dots & 0 & 0 \\ \dots & \dots & \dots & \dots & \dots & \dots & 0 & \dots \\ 0 & 0 & \dots & \cos\theta & 0 & \dots & 0 & \sin\theta \\ 0 & 0 & \dots & 0 & 1 & \dots & 0 & 0 \\ \dots & \dots & \dots & \dots & \dots & \dots & 1 & \dots \\ 0 & 0 & \dots & -\sin\theta & 0 & \dots & 0 & \cos\theta \end{array} \right) $$ with property \( {\bf S^{T}} = {\bf S^{-1}} \). It performs a plane rotation around an angle \( \theta \) in the Euclidean $n-$dimensional space.

It means that its matrix elements that differ from zero are given by $$ s_{kk}= s_{ll}=\cos\theta, s_{kl}=-s_{lk}= -\sin\theta, s_{ii}=1\hspace{0.5cm} i\ne k \hspace{0.5cm} i \ne l, $$ A similarity transformation $$ {\bf B}= {\bf S}^T {\bf A}{\bf S}, $$ results in $$ \begin{eqnarray*} b_{ik} &=& a_{ik}\cos\theta - a_{il}\sin\theta , i \ne k, i \ne l \\ b_{il} &=& a_{il}\cos\theta + a_{ik}\sin\theta , i \ne k, i \ne l \nonumber\\ b_{kk} &=& a_{kk}\cos^2\theta - 2a_{kl}\cos\theta \sin\theta +a_{ll}\sin^2\theta\nonumber\\ b_{ll} &=& a_{ll}\cos^2\theta +2a_{kl}\cos\theta sin\theta +a_{kk}\sin^2\theta\nonumber\\ b_{kl} &=& (a_{kk}-a_{ll})\cos\theta \sin\theta +a_{kl}(\cos^2\theta-\sin^2\theta)\nonumber \end{eqnarray*} $$ The angle \( \theta \) is arbitrary. The recipe is to choose \( \theta \) so that all non-diagonal matrix elements \( b_{kl} \) become zero.

The main idea is thus to reduce systematically the norm of the off-diagonal matrix elements of a matrix \( {\bf A} \) $$ \mathrm{off}({\bf A}) = \sqrt{\sum_{i=1}^n\sum_{j=1,j\ne i}^n a_{ij}^2}. $$ To demonstrate the algorithm, we consider the simple \( 2\times 2 \) similarity transformation of the full matrix. The matrix is symmetric, we single out $ 1\le k < l \le n$ and use the abbreviations \( c=\cos\theta \) and \( s=\sin\theta \) to obtain $$ \left( \begin{array}{cc} b_{kk} & 0 \\ 0 & b_{ll} \\\end{array} \right) = \left( \begin{array}{cc} c & -s \\ s &c \\\end{array} \right) \left( \begin{array}{cc} a_{kk} & a_{kl} \\ a_{lk} &a_{ll} \\\end{array} \right) \left( \begin{array}{cc} c & s \\ -s & c \\\end{array} \right). $$

We require that the non-diagonal matrix elements \( b_{kl}=b_{lk}=0 \), implying that $$ a_{kl}(c^2-s^2)+(a_{kk}-a_{ll})cs = b_{kl} = 0. $$ If \( a_{kl}=0 \) one sees immediately that \( \cos\theta = 1 \) and \( \sin\theta=0 \).

The Frobenius norm of an orthogonal transformation is always preserved. The Frobenius norm is defined as $$ \mathrm{norm}({\bf A})_F = \sqrt{\sum_{i=1}^n\sum_{j=1}^n |a_{ij}|^2}. $$ This means that for our \( 2\times 2 \) case we have $$ 2a_{kl}^2+a_{kk}^2+a_{ll}^2 = b_{kk}^2+b_{ll}^2, $$ which leads to $$ \mathrm{off}({\bf B})^2 = \mathrm{norm}({\bf B})_F^2-\sum_{i=1}^nb_{ii}^2=\mathrm{off}({\bf A})^2-2a_{kl}^2, $$ since $$ \mathrm{norm}({\bf B})_F^2-\sum_{i=1}^nb_{ii}^2=\mathrm{norm}({\bf A})_F^2-\sum_{i=1}^na_{ii}^2+(a_{kk}^2+a_{ll}^2 -b_{kk}^2-b_{ll}^2). $$ This results means that the matrix \( {\bf A} \) moves closer to diagonal form for each transformation.

Defining the quantities \( \tan\theta = t= s/c \) and $$\cot 2\theta=\tau = \frac{a_{ll}-a_{kk}}{2a_{kl}}, $$ we obtain the quadratic equation (using \( \cot 2\theta=1/2(\cot \theta-\tan\theta) \) $$ t^2+2\tau t-1= 0, $$ resulting in $$ t = -\tau \pm \sqrt{1+\tau^2}, $$ and \( c \) and \( s \) are easily obtained via $$ c = \frac{1}{\sqrt{1+t^2}}, $$ and \( s=tc \). Choosing \( t \) to be the smaller of the roots ensures that \( |\theta| \le \pi/4 \) and has the effect of minimizing the difference between the matrices \( {\bf B} \) and \( {\bf A} \) since $$ \mathrm{norm}({\bf B}-{\bf A})_F^2=4(1-c)\sum_{i=1,i\ne k,l}^n(a_{ik}^2+a_{il}^2) +\frac{2a_{kl}^2}{c^2}. $$

• Choose a tolerance \( \epsilon \), making it a small number, typically \( 10^{-8} \) or smaller.
• Setup a while test where one compares the norm of the newly computed off-diagonal

matrix elements \[ \mathrm{off}({\bf A}) = \sqrt{\sum_{i=1}^n\sum_{j=1,j\ne i}^n a_{ij}^2} > \epsilon. \]

• Now choose the matrix elements \( a_{kl} \) so that we have those with largest value, that is \( |a_{kl}|=\mathrm{max}_{i\ne j} |a_{ij}| \).

Compute thereafter \( \tau = (a_{ll}-a_{kk})/2a_{kl} \), \( \tan\theta \), \( \cos\theta \) and \( \sin\theta \).

• Compute thereafter the similarity transformation for this set of values \( (k,l) \), obtaining the new matrix \( {\bf B}= {\bf S}(k,l,\theta)^T {\bf A}{\bf S}(k,l,\theta) \).
• Compute the new norm of the off-diagonal matrix elements and continue till you have satisfied \( \mathrm{off}({\bf B}) \le \epsilon \)

The convergence rate of the Jacobi method is however poor, one needs typically \( 3n^2-5n^2 \) rotations and each rotation requires \( 4n \) operations, resulting in a total of \( 12n^3-20n^3 \) operations in order to zero out non-diagonal matrix elements.

We specialize to a symmetric \( 3\times 3 \) matrix \( {\bf A} \). We start the process as follows (assuming that \( a_{23}=a_{32} \) is the largest non-diagonal) with \( c=\cos{\theta} \) and \( s=\sin{\theta} \) $$ {\bf B} = \left( \begin{array}{ccc} 1 & 0 & 0 \\ 0 & c & -s \\ 0 & s & c \end{array} \right)\left( \begin{array}{ccc} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \end{array} \right) \left( \begin{array}{ccc} 1 & 0 & 0 \\ 0 & c & s \\ 0 & -s & c \end{array} \right). $$ We will choose the angle \( \theta \) in order to have \( a_{23}=a_{32}=0 \). We get (symmetric matrix) $$ {\bf B} =\left( \begin{array}{ccc} a_{11} & a_{12}c -a_{13}s& a_{12}s+a_{13}c \\ a_{12}c -a_{13}s & a_{22}c^2+a_{33}s^2 -2a_{23}sc& (a_{22}-a_{33})sc +a_{23}(c^2-s^2) \\ a_{12}s+a_{13}c & (a_{22}-a_{33})sc +a_{23}(c^2-s^2) & a_{22}s^2+a_{33}c^2 +2a_{23}sc \end{array} \right). $$ Note that \( a_{11} \) is unchanged! As it should.

We have $$ {\bf B} =\left( \begin{array}{ccc} a_{11} & a_{12}c -a_{13}s& a_{12}s+a_{13}c \\ a_{12}c -a_{13}s & a_{22}c^2+a_{33}s^2 -2a_{23}sc& (a_{22}-a_{33})sc +a_{23}(c^2-s^2) \\ a_{12}s+a_{13}c & (a_{22}-a_{33})sc +a_{23}(c^2-s^2) & a_{22}s^2+a_{33}c^2 +2a_{23}sc \end{array} \right). $$ or $$ \begin{eqnarray*} b_{11} &=& a_{11} \\ b_{12} &=& a_{12}\cos\theta - a_{13}\sin\theta , 1 \ne 2, 1 \ne 3 \\ b_{13} &=& a_{13}\cos\theta + a_{12}\sin\theta , 1 \ne 2, 1 \ne 3 \nonumber\\ b_{22} &=& a_{22}\cos^2\theta - 2a_{23}\cos\theta \sin\theta +a_{33}\sin^2\theta\nonumber\\ b_{33} &=& a_{33}\cos^2\theta +2a_{23}\cos\theta \sin\theta +a_{22}\sin^2\theta\nonumber\\ b_{23} &=& (a_{22}-a_{33})\cos\theta \sin\theta +a_{23}(\cos^2\theta-\sin^2\theta)\nonumber \end{eqnarray*} $$ We will fix the angle \( \theta \) so that \( b_{23}=0 \).

We get then a new matrix $$ {\bf B} =\left( \begin{array}{ccc} b_{11} & b_{12}& b_{13} \\ b_{12}& b_{22}& 0 \\ b_{13}& 0& a_{33} \end{array} \right). $$ We repeat then assuming that \( b_{12} \) is the largest non-diagonal matrix element and get a new matrix $$ {\bf C} = \left( \begin{array}{ccc} c & -s & 0 \\ s & c & 0 \\ 0 & 0 & 1 \end{array} \right)\left( \begin{array}{ccc} b_{11} & b_{12} & b_{13} \\ b_{12} & b_{22} & 0 \\ b_{13} & 0 & b_{33} \end{array} \right) \left( \begin{array}{ccc} c & s & 0 \\ -s & c & 0 \\ 0 & 0 & 1 \end{array} \right). $$ We continue this process till all non-diagonal matrix elements are zero (ideally). You will notice that performing the above operations that the matrix element \( b_{23} \) which was previous zero becomes different from zero. This is one of the problems which slows down the jacobi procedure.

The more general expression for the new matrix elements are $$ \begin{eqnarray*} b_{ii} &=& a_{ii}, i \ne k, i \ne l \\ b_{ik} &=& a_{ik}\cos\theta - a_{il}\sin\theta , i \ne k, i \ne l \\ b_{il} &=& a_{il}\cos\theta + a_{ik}\sin\theta , i \ne k, i \ne l \nonumber\\ b_{kk} &=& a_{kk}\cos^2\theta - 2a_{kl}\cos\theta \sin\theta +a_{ll}\sin^2\theta\nonumber\\ b_{ll} &=& a_{ll}\cos^2\theta +2a_{kl}\cos\theta \sin\theta +a_{kk}\sin^2\theta\nonumber\\ b_{kl} &=& (a_{kk}-a_{ll})\cos\theta \sin\theta +a_{kl}(\cos^2\theta-\sin^2\theta)\nonumber \end{eqnarray*} $$ This is what we will need to code.

Code example.

//  we have defined a matrix A and a matrix R for the eigenvector, both of dim n x n
//  The final matrix R has the eigenvectors in its row elements, it is set to one
//  for the diagonal elements in the beginning, zero else.
....
double tolerance = 1.0E-10; 
int iterations = 0;
while ( maxnondiag > tolerance && iterations <= maxiter)
{
   int p, q;
   maxnondiag  = offdiag(A, p, q, n);
   Jacobi_rotate(A, R, p, q, n);
   iterations++;
}
...

Finding the max nondiagonal element

//  the offdiag function, using Armadillo
double offdiag(mat A, int p, int q, int n);
{
   double max;
   for (int i = 0; i < n; ++i)
   {
       for ( int j = i+1; j < n; ++j)
       {
           double aij = fabs(A(i,j));
           if ( aij > max)
           { 
              max = aij;  p = i; q = j;
           }
       }
   }
   return max;
}
// more statements

Finding the new matrix elements

void Jacobi_rotate ( mat A, mat R, int k, int l, int n )
{
  double s, c;
  if ( A(k,l) != 0.0 ) {
    double t, tau;
    tau = (A(l,l) - A(k,k))/(2*A(k,l));
    
    if ( tau >= 0 ) {
      t = 1.0/(tau + sqrt(1.0 + tau*tau));
    } else {
      t = -1.0/(-tau +sqrt(1.0 + tau*tau));
    }
    
    c = 1/sqrt(1+t*t);
    s = c*t;
  } else {
    c = 1.0;
    s = 0.0;
  }
  double a_kk, a_ll, a_ik, a_il, r_ik, r_il;
  a_kk = A(k,k);
  a_ll = A(l,l);
  A(k,k) = c*c*a_kk - 2.0*c*s*A(k,l) + s*s*a_ll;
  A(l,l) = s*s*a_kk + 2.0*c*s*A(k,l) + c*c*a_ll;
  A(k,l) = 0.0;  // hard-coding non-diagonal elements by hand
  A(l,k) = 0.0;  // same here
  for ( int i = 0; i < n; i++ ) {
    if ( i != k && i != l ) {
      a_ik = A(i,k);
      a_il = A(i,l);
      A(i,k) = c*a_ik - s*a_il;
      A(k,i) = A(i,k);
      A(i,l) = c*a_il + s*a_ik;
      A(l,i) = A(i,l);
    }
//  And finally the new eigenvectors
    r_ik = R(i,k);
    r_il = R(i,l);

    R(i,k) = c*r_ik - s*r_il;
    R(i,l) = c*r_il + s*r_ik;
  }
  return;
} // end of function jacobi_rotate

The drawbacks with Jacobi's method are rather obvious, with perhaps the most negative feature being the fact that we cannot tell * a priori* how many transformations are needed. Can we do better? The answer to this is yes and is given by a clever algorithm outlined by Householder. It was ranked among the top ten algorithms in the previous century. We will discuss this algorithm in more detail during week 38.

The first step consists in finding an orthogonal matrix \( {\bf S} \) which is the product of \( (n-2) \) orthogonal matrices $$ {\bf S}={\bf S}_1{\bf S}_2\dots{\bf S}_{n-2}, $$ each of which successively transforms one row and one column of \( {\bf A} \) into the required tridiagonal form. Only \( n-2 \) transformations are required, since the last two elements are already in tridiagonal form. In order to determine each \( {\bf S_i} \) let us see what happens after the first multiplication, namely, $$ {\bf S}_1^T{\bf A}{\bf S}_1= \left( \begin{array}{ccccccc} a_{11} & e_1 & 0 & 0 & \dots &0 & 0 \\ e_1 & a'_{22} &a'_{23} & \dots & \dots &\dots &a'_{2n} \\ 0 & a'_{32} &a'_{33} & \dots & \dots &\dots &a'_{3n} \\ 0 & \dots &\dots & \dots & \dots &\dots & \\ 0 & a'_{n2} &a'_{n3} & \dots & \dots &\dots &a'_{nn} \\ \end{array} \right) $$ where the primed quantities represent a matrix \( {\bf A}' \) of dimension \( n-1 \) which will subsequentely be transformed by \( {\bf S_2} \).

The factor \( e_1 \) is a possibly non-vanishing element. The next transformation produced by \( {\bf S_2} \) has the same effect as \( {\bf S_1} \) but now on the submatirx \( {\bf A^{'}} \) only $$ \left ({\bf S}_{1}{\bf S}_{2} \right )^{T} {\bf A}{\bf S}_{1} {\bf S}_{2} = \left( \begin{array}{ccccccc} a_{11} & e_1 & 0 & 0 & \dots &0 & 0 \\ e_1 & a'_{22} &e_2 & 0 & \dots &\dots &0 \\ 0 & e_2 &a''_{33} & \dots & \dots &\dots &a''_{3n} \\ 0 & \dots &\dots & \dots & \dots &\dots & \\ 0 & 0 &a''_{n3} & \dots & \dots &\dots &a''_{nn} \\ \end{array} \right) $$ Note that the effective size of the matrix on which we apply the transformation reduces for every new step. In the previous Jacobi method each similarity transformation is in principle performed on the full size of the original matrix.

After a series of such transformations, we end with a set of diagonal matrix elements $$ a_{11}, a'_{22}, a''_{33}\dots a^{n-1}_{nn}, $$ and off-diagonal matrix elements $$ e_1, e_2,e_3, \dots, e_{n-1}. $$ The resulting matrix reads $$ {\bf S}^{T} {\bf A} {\bf S} = \left( \begin{array}{ccccccc} a_{11} & e_1 & 0 & 0 & \dots &0 & 0 \\ e_1 & a'_{22} & e_2 & 0 & \dots &0 &0 \\ 0 & e_2 & a''_{33} & e_3 &0 &\dots & 0\\ \dots & \dots & \dots & \dots &\dots &\dots & \dots\\ 0 & \dots & \dots & \dots &\dots &a^{(n-1)}_{n-2} & e_{n-1}\\ 0 & \dots & \dots & \dots &\dots &e_{n-1} & a^{(n-1)}_{n-1} \end{array} \right) . $$

It remains to find a recipe for determining the transformation \( {\bf S}_n \). We illustrate the method for \( {\bf S}_1 \) which we assume takes the form $$ {\bf S_{1}} = \left( \begin{array}{cc} 1 & {\bf 0^{T}} \\ {\bf 0}& {\bf P} \end{array} \right), $$ with \( {\bf 0^{T}} \) being a zero row vector, \( {\bf 0^{T}} = \{0,0,\cdots\} \) of dimension \( (n-1) \). The matrix \( {\bf P} \) is symmetric with dimension (\( (n-1) \times (n-1) \)) satisfying \( {\bf P}^2={\bf I} \) and \( {\bf P}^T={\bf P} \). A possible choice which fullfils the latter two requirements is $$ {\bf P}={\bf I}-2{\bf u}{\bf u}^T, $$ where \( {\bf I} \) is the \( (n-1) \) unity matrix and \( {\bf u} \) is an \( n-1 \) column vector with norm \( {\bf u}^T{\bf u} \) (inner product).

Note that \( {\bf u}{\bf u}^T \) is an outer product giving a matrix of dimension (\( (n-1) \times (n-1) \)). Each matrix element of \( {\bf P} \) then reads $$ P_{ij}=\delta_{ij}-2u_iu_j, $$ where \( i \) and \( j \) range from \( 1 \) to \( n-1 \). Applying the transformation \( {\bf S}_1 \) results in $$ {\bf S}_1^T{\bf A}{\bf S}_1 = \left( \begin{array}{cc} a_{11} & ({\bf Pv})^T \\ {\bf Pv}& {\bf A}' \end{array} \right) , $$ where \( {\bf v^{T}} = \{a_{21}, a_{31},\cdots, a_{n1}\} \) and {\bf P} must satisfy (\( {\bf Pv})^{T} = \{k, 0, 0,\cdots \} \). Then $$ \begin{equation} {\bf Pv} = {\bf v} -2{\bf u}( {\bf u}^T{\bf v})= k {\bf e}, \label{eq:palpha} \end{equation} $$ with \( {\bf e^{T}} = \{1,0,0,\dots 0\} \).

Solving the latter equation gives us \( {\bf u} \) and thus the needed transformation \( {\bf P} \). We do first however need to compute the scalar \( k \) by taking the scalar product of the last equation with its transpose and using the fact that \( {\bf P}^2={\bf I} \). We get then $$ ({\bf Pv})^T{\bf Pv} = k^{2} = {\bf v}^T{\bf v}= |v|^2 = \sum_{i=2}^{n}a_{i1}^2, $$ which determines the constant $ k = \pm v$.

Now we can rewrite Eq. \eqref{eq:palpha} as $$ {\bf v} - k{\bf e} = 2{\bf u}( {\bf u}^T{\bf v}), $$ and taking the scalar product of this equation with itself and obtain $$ \begin{equation} 2( {\bf u}^T{\bf v})^2=(v^2\pm a_{21}v), \label{eq:pmalpha} \end{equation} $$ which finally determines $$ {\bf u}=\frac{{\bf v}-k{\bf e}}{2( {\bf u}^T{\bf v})}. $$ In solving Eq. \eqref{eq:pmalpha} great care has to be exercised so as to choose those values which make the right-hand largest in order to avoid loss of numerical precision. The above steps are then repeated for every transformations till we have a tridiagonal matrix suitable for obtaining the eigenvalues.

Our Householder transformation has given us a tridiagonal matrix. We discuss here how one can use Householder's iterative procedure to obtain the eigenvalues. Let us specialize to a \( 4\times 4 \) matrix. The tridiagonal matrix takes the form $$ {\bf A} = \left( \begin{array}{cccc} d_{1} & e_{1} & 0 & 0 \\ e_{1} & d_{2} & e_{2} & 0 \\ 0 & e_{2} & d_{3} & e_{3} \\ 0 & 0 & e_{3} & d_{4} \end{array} \right). $$ As a first observation, if any of the elements \( e_{i} \) are zero the matrix can be separated into smaller pieces before diagonalization. Specifically, if \( e_{1} = 0 \) then \( d_{1} \) is an eigenvalue.

Thus, let us introduce a transformation \( {\bf S_{1}} \) which operates like $$ {\bf S_{1}} = \left( \begin{array}{cccc} \cos \theta & 0 & 0 & \sin \theta\\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ \cos \theta & 0 & 0 & \cos \theta \end{array} \right) $$

Then the similarity transformation $$ {\bf S_{1}^{T} A S_{1}} = {\bf A'} = \left( \begin{array}{cccc} d'_{1} & e'_{1} & 0 & 0 \\ e'_{1} & d_{2} & e_{2} & 0 \\ 0 & e_{2} & d_{3} & e'{3} \\ 0 & 0 & e'_{3} & d'_{4} \end{array} \right) $$ produces a matrix where the primed elements in \( {\bf A'} \) have been changed by the transformation whereas the unprimed elements are unchanged. If we now choose \( \theta \) to give the element \( a_{21}^{'} = e^{'}= 0 \) then we have the first eigenvalue \( = a_{11}^{'} = d_{1}^{'} \). (This is actually what you are doing in project 2!!)

This procedure can be continued on the remaining three-dimensional submatrix for the next eigenvalue. Thus after few transformations we have the wanted diagonal form.

What we see here is just a special case of the more general procedure developed by Francis in two articles in 1961 and 1962.

The algorithm is based on the so-called QR method (or just QR-algorithm). It follows from a theorem by Schur which states that any square matrix can be written out in terms of an orthogonal matrix \( {\bf Q} \) and an upper triangular matrix \( {\bf U} \). Historically R was used instead of U since the wording right triangular matrix was first used. The method is based on an iterative procedure similar to Jacobi's method, by a succession of planar rotations. For a tridiagonal matrix it is simple to carry out in principle, but complicated in detail! We will discuss this in more detail during week 38.

Eigenvalue problems.

• Monday: Brief repetition from last week, with discussion of project 2.
• Mathematical aspects of Jacobi's algorithm
• Discussion of Householder's algorithm
• Tuesday: More on Householder's algorithm and diagonalization of tridiagonal matrices (Givens and Francis' algorithms)
• Discussion of Lanczos' method (also optional part of project 2)
• Lab: discussion of unit tests and how to write a scientific report

Reading assignment this week: chapters 7.1-7.4 of lecture notes

Our Householder transformation has given us a tridiagonal matrix. We discuss here how one can use Jacobi's iterative procedure to obtain the eigenvalues, although it may not be the best approach. Let us specialize to a \( 4\times 4 \) matrix. The tridiagonal matrix takes the form $$ {\bf A} = \left( \begin{array}{cccc} d_{1} & e_{1} & 0 & 0 \\ e_{1} & d_{2} & e_{2} & 0 \\ 0 & e_{2} & d_{3} & e_{3} \\ 0 & 0 & e_{3} & d_{4} \end{array} \right). $$ As a first observation, if any of the elements \( e_{i} \) are zero the matrix can be separated into smaller pieces before diagonalization. Specifically, if \( e_{1} = 0 \) then \( d_{1} \) is an eigenvalue.

Thus, let us introduce a transformation \( {\bf S_{1}} \) which operates like $$ {\bf S_{1}} = \left( \begin{array}{cccc} \cos \theta & 0 & 0 & \sin \theta\\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ \cos \theta & 0 & 0 & \cos \theta \end{array} \right) $$ Then the similarity transformation $$ {\bf S_{1}^{T} A S_{1}} = {\bf A'} = \left( \begin{array}{cccc} d'_{1} & e'_{1} & 0 & 0 \\ e'_{1} & d_{2} & e_{2} & 0 \\ 0 & e_{2} & d_{3} & e'{3} \\ 0 & 0 & e'_{3} & d'_{4} \end{array} \right) $$ produces a matrix where the primed elements in \( {\bf A'} \) have been changed by the transformation whereas the unprimed elements are unchanged. If we now choose \( \theta \) to give the element \( a_{21}^{'} = e^{'}= 0 \) then we have the first eigenvalue \( = a_{11}^{'} = d_{1}^{'} \).

This procedure can be continued on the remaining three-dimensional submatrix for the next eigenvalue. Thus after few transformations we have the wanted diagonal form.

What we see here is just a special case of the more general procedure developed by Francis in two articles in 1961 and 1962. Using Jacobi's method is not very efficient ether.

The algorithm is based on the so-called {\bf QR} method (or just {\bf QR}-algorithm). It follows from a theorem by Schur which states that any square matrix can be written out in terms of an orthogonal matrix \( \hat{Q} \) and an upper triangular matrix \( \hat{U} \). Historically \( R \) was used instead of \( U \) since the wording right triangular matrix was first used.

The method is based on an iterative procedure similar to Jacobi's method, by a succession of planar rotations. For a tridiagonal matrix it is simple to carry out in principle, but complicated in detail!

Schur's theorem $$ \hat{A} = \hat{Q}\hat{U}, $$ is used to rewrite any square matrix into a unitary matrix times an upper triangular matrix. We say that a square matrix is similar to a triangular matrix.

Householder's algorithm which we have derived is just a special case of the general Householder algorithm. For a symmetric square matrix we obtain a tridiagonal matrix.

There is a corollary to Schur's theorem which states that every Hermitian matrix is unitarily similar to a diagonal matrix.

It follows that we can define a new matrix $$ \hat{A}\hat{Q} = \hat{Q}\hat{U}\hat{Q}, $$ and multiply from the left with \( \hat{Q}^{-1} \) we get $$ \hat{Q}^{-1}\hat{A}\hat{Q} = \hat{B}=\hat{U}\hat{Q}, $$ where the matrix \( \hat{B} \) is a similarity transformation of \( \hat{A} \) and has the same eigenvalues as \( \hat{B} \).

Suppose \( \hat{A} \) is the triangular matrix we obtained after the Householder transformation, $$ \hat{A} = \hat{Q}\hat{U}, $$ and multiply from the left with \( \hat{Q}^{-1} \) resulting in $$ \hat{Q}^{-1}\hat{A} = \hat{U}. $$ Suppose that \( \hat{Q} \) consists of a series of planar Jacobi like rotations acting on sub blocks of \( \hat{A} \) so that all elements below the diagonal are zeroed out $$ \hat{Q}=\hat{R}_{12}\hat{R}_{23}\dots\hat{R}_{n-1,n}. $$

A transformation of the type \( \hat{R}_{12} \) looks like $$ \hat{R}_{12} = \left( \begin{array}{ccccccccc} c&s &0 &0 &0 & \dots &0 & 0 & 0\\ -s&c &0 &0 &0 & \dots &0 & 0 & 0\\ 0&0 &1 &0 &0 & \dots &0 & 0 & 0\\ \dots&\dots &\dots &\dots &\dots &\dots \\ 0&0 &0 & 0 & 0 & \dots &1 &0 &0 \\ 0&0 &0 & 0 & 0 & \dots &0 &1 &0 \\ 0&0 &0 & 0 & 0 & \dots &0 &0 & 1 \end{array} \right) $$

The matrix \( \hat{U} \) takes then the form $$ \hat{U} = \left( \begin{array}{ccccccccc} x&x &x &0 &0 & \dots &0 & 0 & 0\\ 0&x &x &x &0 & \dots &0 & 0 & 0\\ 0&0 &x &x &x & \dots &0 & 0 & 0\\ \dots&\dots &\dots &\dots &\dots &\dots \\ 0&0 &0 & 0 & 0 & \dots &x &x &x \\ 0&0 &0 & 0 & 0 & \dots &0 &x &x \\ 0&0 &0 & 0 & 0 & \dots &0 &0 & x \end{array} \right) $$ which has a second superdiagonal.

We have now found \( \hat{Q} \) and \( \hat{U} \) and this allows us to find the matrix \( \hat{B} \) which is, due to Schur's theorem, unitarily similar to a triangular matrix (upper in our case) since we have that $$ \hat{Q}^{-1}\hat{A}\hat{Q} = \hat{B}, $$ from Schur's theorem the matrix \( \hat{B} \) is triangular and the eigenvalues the same as those of \( \hat{A} \) and are given by the diagonal matrix elements of \( \hat{B} \). Why?

Our matrix \( \hat{B}=\hat{U}\hat{Q} \).

The matrix \( \hat{A} \) is transformed into a tridiagonal form and the last step is to transform it into a diagonal matrix giving the eigenvalues on the diagonal.

The eigenvalues of a matrix can be obtained using the characteristic polynomial $$ P(\lambda) = det(\lambda{\bf I}-{\bf A})= \prod_{i=1}^{n}\left(\lambda_i-\lambda\right), $$ which rewritten in matrix form reads $$ P(\lambda)= \left( \begin{array}{ccccccc} d_1-\lambda & e_1 & 0 & 0 & \dots &0 & 0 \\ e_1 & d_2-\lambda & e_2 & 0 & \dots &0 &0 \\ 0 & e_2 & d_3-\lambda & e_3 &0 &\dots & 0\\ \dots & \dots & \dots & \dots &\dots &\dots & \dots\\ 0 & \dots & \dots & \dots &\dots &d_{N_{\mathrm{step}}-2}-\lambda & e_{N_{\mathrm{step}}-1}\\ 0 & \dots & \dots & \dots &\dots &e_{N_{\mathrm{step}}-1} & d_{N_{\mathrm{step}}-1}-\lambda \end{array} \right) $$ We can solve this equation in an iterative manner. We let \( P_k(\lambda) \) be the value of \( k \) subdeterminant of the above matrix of dimension \( n\times n \). The polynomial \( P_k(\lambda) \) is clearly a polynomial of degree \( k \). Starting with \( P_1(\lambda) \) we have \( P_1(\lambda)=d_1-\lambda \). The next polynomial reads \( P_2(\lambda)=(d_2-\lambda)P_1(\lambda)-e_1^2 \). By expanding the determinant for \( P_k(\lambda) \) in terms of the minors of the $n$th column we arrive at the recursion relation \[ P_k(\lambda)=(d_k-\lambda)P_{k-1}(\lambda)-e_{k-1}^2P_{k-2}(\lambda). \] Together with the starting values \( P_1(\lambda) \) and \( P_2(\lambda) \) and good root searching methods we arrive at an efficient computational scheme for finding the roots of \( P_n(\lambda) \). However, for large matrices this algorithm is rather inefficient and time-consuming.

Basic features with a real symmetric matrix (and normally huge \( n> 10^6 \) and sparse) \( \hat{A} \) of dimension \( n\times n \):

• Lanczos' algorithm generates a sequence of real tridiagonal matrices \( T_k \) of dimension \( k\times k \) with \( k\le n \), with the property that the extremal eigenvalues of \( T_k \) are progressively better estimates of \( \hat{A} \)' extremal eigenvalues.* The method converges to the extremal eigenvalues.
• The similarity transformation is

$$ \hat{T}= \hat{Q}^{T}\hat{A}\hat{Q}, $$ with the first vector \( \hat{Q}\hat{e}_1=\hat{q}_1 \).

We are going to solve iteratively $$ \hat{T}= \hat{Q}^{T}\hat{A}\hat{Q}, $$ with the first vector \( \hat{Q}\hat{e}_1=\hat{q}_1 \). We can write out the matrix \( \hat{Q} \) in terms of its column vectors $$ \hat{Q}=\left[\hat{q}_1\hat{q}_2\dots\hat{q}_n\right]. $$

The matrix $$ \hat{T}= \hat{Q}^{T}\hat{A}\hat{Q}, $$ can be written as $$ \hat{T} = \left(\begin{array}{cccccc} \alpha_1& \beta_1 & 0 &\dots & \dots &0 \\ \beta_1 & \alpha_2 & \beta_2 &0 &\dots &0 \\ 0& \beta_2 & \alpha_3 & \beta_3 & \dots &0 \\ \dots& \dots & \dots &\dots &\dots & 0 \\ \dots& & &\beta_{n-2} &\alpha_{n-1}& \beta_{n-1} \\ 0& \dots &\dots &0 &\beta_{n-1} & \alpha_{n} \\ \end{array} \right) $$

Using the fact that \( \hat{Q}\hat{Q}^T=\hat{I} \), we can rewrite $$ \hat{T}= \hat{Q}^{T}\hat{A}\hat{Q}, $$ as $$ \hat{Q}\hat{T}= \hat{A}\hat{Q}, $$ and if we equate columns $$ \hat{T} = \left(\begin{array}{cccccc} \alpha_1& \beta_1 & 0 &\dots & \dots &0 \\ \beta_1 & \alpha_2 & \beta_2 &0 &\dots &0 \\ 0& \beta_2 & \alpha_3 & \beta_3 & \dots &0 \\ \dots& \dots & \dots &\dots &\dots & 0 \\ \dots& & &\beta_{n-2} &\alpha_{n-1}& \beta_{n-1} \\ 0& \dots &\dots &0 &\beta_{n-1} & \alpha_{n} \\ \end{array} \right) $$ we obtain $$ \hat{A}\hat{q}_k=\beta_{k-1}\hat{q}_{k-1}+\alpha_k\hat{q}_k+\beta_k\hat{q}_{k+1}. $$

We have thus $$ \hat{A}\hat{q}_k=\beta_{k-1}\hat{q}_{k-1}+\alpha_k\hat{q}_k+\beta_k\hat{q}_{k+1}, $$ with \( \beta_0\hat{q}_0=0 \) for \( k=1:n-1 \). Remember that the vectors \( \hat{q}_k \) are orthornormal and this implies $$ \alpha_k=\hat{q}_k^T\hat{A}\hat{q}_k, $$ and these vectors are called Lanczos vectors.

We have thus $$ \hat{A}\hat{q}_k=\beta_{k-1}\hat{q}_{k-1}+\alpha_k\hat{q}_k+\beta_k\hat{q}_{k+1}, $$ with \( \beta_0\hat{q}_0=0 \) for \( k=1:n-1 \) and $$ \alpha_k=\hat{q}_k^T\hat{A}\hat{q}_k. $$ If $$ \hat{r}_k=(\hat{A}-\alpha_k\hat{I})\hat{q}_k-\beta_{k-1}\hat{q}_{k-1}, $$ is non-zero, then $$ \hat{q}_{k+1}=\hat{r}_{k}/\beta_k, $$ with \( \beta_k=\pm ||\hat{r}_{k}||_2 \).

What should it contain? A typical structure.

• An introduction where you explain the aims and rationale for the physics case and what you have done. At the end of the introduction you should give a brief summary of the structure of the report
• Theoretical models and technicalities. This is the methods section.
• Results and discussion
• Conclusions and perspectives
• Appendix with extra material
• Bibliography

Keep always a good log of what you do.

What should I focus on? Introduction.

You don't need to answer all questions in a chronological order. When you write the introduction you could focus on the following aspects

• Motivate the reader, the first part of the introduction gives always a motivation and tries to give the overarching ideas
• What I have done
• The structure of the report, how it is organized etc

What should I focus on? Methods sections.

• Describe the methods and algorithms
• You need to explain how you implemented the methods and also say something about the structure of your algorithm and present some parts of your code
• You should plug in some calculations to demonstrate your code, such as selected runs used to validate and verify your results. The latter is extremely important!! A reader needs to understand that your code reproduces selected benchmarks and reproduces previous results, either numerical and/or well-known closed form expressions.

What should I focus on? Results.

• Present your results
• Give a critical discussion of your work and place it in the correct context.
• Relate your work to other calculations/studies
• An eventual reader should be able to reproduce your calculations if she/he wants to do so. All input variables should be properly explained.
• Make sure that figures and tables should contain enough information in their captions, axis labels etc so that an eventual reader can gain a first impression of your work by studying figures and tables only.

What should I focus on? Conclusions.

• State your main findings and interpretations
• Try as far as possible to present perspectives for future work
• Try to discuss the pros and cons of the methods and possible improvements

What should I focus on? additional material.

• Additional calculations used to validate the codes
• Selected calculations, these can be listed with few comments
• Listing of the code if you feel this is necessary

You can consider moving parts of the material from the methods section to the appendix. You can also place additional material on your webpage.

What should I focus on? References.

• Give always references to material you base your work on, either scientific articles/reports or books.
• Refer to articles as: name(s) of author(s), journal, volume (boldfaced), page and year in parenthesis.
• Refer to books as: name(s) of author(s), title of book, publisher, place and year, eventual page numbers

Unit Testing is the practice of testing the smallest testable parts, called units, of an application individually and independently to determine if they behave exactly as expected. Unit tests (short code fragments) are usually written such that they can be preformed at any time during the development to continually verify the behavior of the code. In this way, possible bugs will be identified early in the development cycle, making the debugging at later stage much easier. There are many benefits associated with Unit Testing, such as

• It increases confidence in changing and maintaining code. Big changes can be made to the code quickly, since the tests will ensure that everything still is working properly.
• Since the code needs to be modular to make Unit Testing possible, the code will be easier to reuse. This improves the code design.
• Debugging is easier, since when a test fails, only the latest changes need to be debugged.
• Different parts of a project can be tested without the need to wait for the other parts to be available.
• A unit test can serve as a documentation on the functionality of a unit of the code.

Look up the guide on how to install unit tests for c++ at course webpage. This is the version with classes.

#include <unittest++/UnitTest++.h>

class MyMultiplyClass{
public:
    double multiply(double x, double y) {
        return x * y;
    }
};

TEST(MyMath) {
    MyMultiplyClass my;
    CHECK_EQUAL(56, my.multiply(7,8));
}

int main()
{
    return UnitTest::RunAllTests();
}

And without classes

#include <unittest++/UnitTest++.h>


double multiply(double x, double y) {
    return x * y;
}

TEST(MyMath) {
    CHECK_EQUAL(56, multiply(7,8));
}

int main()
{
    return UnitTest::RunAllTests();
} 

For Fortran users, the link at http://sourceforge.net/projects/fortranxunit/ contains a similar software for unit testing.

• Direct solvers such as Gauss elimination and LU decomposition discussed in connection with project 1.
• Iterative solvers such as Basic iterative solvers, Jacobi, Gauss-Seidel, Successive over-relaxation. These methods are easy to parallelize, as we will se later. Much used in solutions of partial differential equations.
• Other iterative methods such as Krylov subspace methods with Generalized minimum residual (GMRES) and Conjugate gradient etc will not be discussed.

It is a simple method for solving $$ \hat{A}{\bf x}={\bf b}, $$ where \( \hat{A} \) is a matrix and \( {\bf x} \) and \( {\bf b} \) are vectors. The vector \( {\bf x} \) is the unknown.

It is an iterative scheme where we start with a guess for the unknown, and after \( k+1 \) iterations we have $$ {\bf x}^{(k+1)}= \hat{D}^{-1}({\bf b}-(\hat{L}+\hat{U}){\bf x}^{(k)}), $$ with \( \hat{A}=\hat{D}+\hat{U}+\hat{L} \) and \( \hat{D} \) being a diagonal matrix, \( \hat{U} \) an upper triangular matrix and \( \hat{L} \) a lower triangular matrix.

If the matrix \( \hat{A} \) is positive definite or diagonally dominant, one can show that this method will always converge to the exact solution.

We can demonstrate Jacobi's method by this \( 4\times 4 \) matrix problem. We assume a guess for the vector elements \( x_i^{(0)} \), a guess which represents our first iteration. The new values are obtained by substitution $$ \begin{eqnarray} x_1^{(1)} =&(b_1-a_{12}x_2^{(0)} -a_{13}x_3^{(0)} - a_{14}x_4^{(0)})/a_{11} \nonumber \\ x_2^{(1)} =&(b_2-a_{21}x_1^{(0)} - a_{23}x_3^{(0)} - a_{24}x_4^{(0)})/a_{22} \nonumber \\ x_3^{(1)} =&(b_3- a_{31}x_1^{(0)} -a_{32}x_2^{(0)} -a_{34}x_4^{(0)})/a_{33} \nonumber \\ x_4^{(1)}=&(b_4-a_{41}x_1^{(0)} -a_{42}x_2^{(0)} - a_{43}x_3^{(0)})/a_{44}, \nonumber \end{eqnarray} $$ which after \( k+1 \) iterations reads $$ \begin{eqnarray} x_1^{(k+1)} =&(b_1-a_{12}x_2^{(k)} -a_{13}x_3^{(k)} - a_{14}x_4^{(k)})/a_{11} \nonumber \\ x_2^{(k+1)} =&(b_2-a_{21}x_1^{(k)} - a_{23}x_3^{(k)} - a_{24}x_4^{(k)})/a_{22} \nonumber \\ x_3^{(k+1)} =&(b_3- a_{31}x_1^{(k)} -a_{32}x_2^{(k)} -a_{34}x_4^{(k)})/a_{33} \nonumber \\ x_4^{(k+1)}=&(b_4-a_{41}x_1^{(k)} -a_{42}x_2^{(k)} - a_{43}x_3^{(k)})/a_{44}, \nonumber \end{eqnarray} $$

We can generalize the above equations to $$ x_i^{(k+1)}=(b_i-\sum_{j=1, j\ne i}^{n}a_{ij}x_j^{(k)})/a_{ii} $$ or in an even more compact form as $$ {\bf x}^{(k+1)}= \hat{D}^{-1}({\bf b}-(\hat{L}+\hat{U}){\bf x}^{(k)}), $$ with \( \hat{A}=\hat{D}+\hat{U}+\hat{L} \) and \( \hat{D} \) being a diagonal matrix, \( \hat{U} \) an upper triangular matrix and \( \hat{L} \) a lower triangular matrix.

Our \( 4\times 4 \) matrix problem $$ \begin{eqnarray} x_1^{(k+1)} =&(b_1-a_{12}x_2^{(k)} -a_{13}x_3^{(k)} - a_{14}x_4^{(k)})/a_{11} \nonumber \\ x_2^{(k+1)} =&(b_2-a_{21}x_1^{(k)} - a_{23}x_3^{(k)} - a_{24}x_4^{(k)})/a_{22} \nonumber \\ x_3^{(k+1)} =&(b_3- a_{31}x_1^{(k)} -a_{32}x_2^{(k)} -a_{34}x_4^{(k)})/a_{33} \nonumber \\ x_4^{(k+1)}=&(b_4-a_{41}x_1^{(k)} -a_{42}x_2^{(k)} - a_{43}x_3^{(k)})/a_{44}, \nonumber \end{eqnarray} $$ can be rewritten as $$ \begin{eqnarray} x_1^{(k+1)} =&(b_1-a_{12}x_2^{(k)} -a_{13}x_3^{(k)} - a_{14}x_4^{(k)})/a_{11} \nonumber \\ x_2^{(k+1)} =&(b_2-a_{21}x_1^{(k+1)} - a_{23}x_3^{(k)} - a_{24}x_4^{(k)})/a_{22} \nonumber \\ x_3^{(k+1)} =&(b_3- a_{31}x_1^{(k+1)} -a_{32}x_2^{(k+1)} -a_{34}x_4^{(k)})/a_{33} \nonumber \\ x_4^{(k+1)}=&(b_4-a_{41}x_1^{(k+1)} -a_{42}x_2^{(k+1)} - a_{43}x_3^{(k+1)})/a_{44}, \nonumber \end{eqnarray} $$ which allows us to utilize the preceding solution (forward substitution). This improves normally the convergence behavior and leads to the Gauss-Seidel method!

We can generalize $$ \begin{eqnarray} x_1^{(k+1)} =&(b_1-a_{12}x_2^{(k)} -a_{13}x_3^{(k)} - a_{14}x_4^{(k)})/a_{11} \nonumber \\ x_2^{(k+1)} =&(b_2-a_{21}x_1^{(k+1)} - a_{23}x_3^{(k)} - a_{24}x_4^{(k)})/a_{22} \nonumber \\ x_3^{(k+1)} =&(b_3- a_{31}x_1^{(k+1)} -a_{32}x_2^{(k+1)} -a_{34}x_4^{(k)})/a_{33} \nonumber \\ x_4^{(k+1)}=&(b_4-a_{41}x_1^{(k+1)} -a_{42}x_2^{(k+1)} - a_{43}x_3^{(k+1)})/a_{44}, \nonumber \end{eqnarray} $$ to the following form $$ x^{(k+1)}_i = \frac{1}{a_{ii}} \left(b_i - \sum_{j > i}a_{ij}x^{(k)}_j - \sum_{j < i}a_{ij}x^{(k+1)}_j \right),\quad i=1,2,\ldots,n. $$ The procedure is generally continued until the changes made by an iteration are below some tolerance.

The convergence properties of the Jacobi method and the Gauss-Seidel method are dependent on the matrix \( \hat{A} \). These methods converge when the matrix is symmetric positive-definite, or is strictly or irreducibly diagonally dominant. Both methods sometimes converge even if these conditions are not satisfied.

Given a square system of n linear equations with unknown \( \mathbf x \): $$ \hat{A}\mathbf x = \mathbf b $$ where $$ \hat{A}=\begin{bmatrix} a_{11} & a_{12} & \cdots & a_{1n} \\ a_{21} & a_{22} & \cdots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\a_{n1} & a_{n2} & \cdots & a_{nn} \end{bmatrix}, \qquad \mathbf{x} = \begin{bmatrix} x_{1} \\ x_2 \\ \vdots \\ x_n \end{bmatrix} , \qquad \mathbf{b} = \begin{bmatrix} b_{1} \\ b_2 \\ \vdots \\ b_n \end{bmatrix}. $$

Then A can be decomposed into a diagonal component D, and strictly lower and upper triangular components L and U: $$ \hat{A} =\hat{D} + \hat{L} + \hat{U}, $$ where $$ D = \begin{bmatrix} a_{11} & 0 & \cdots & 0 \\ 0 & a_{22} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\0 & 0 & \cdots & a_{nn} \end{bmatrix}, \quad L = \begin{bmatrix} 0 & 0 & \cdots & 0 \\ a_{21} & 0 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\a_{n1} & a_{n2} & \cdots & 0 \end{bmatrix}, \quad U = \begin{bmatrix} 0 & a_{12} & \cdots & a_{1n} \\ 0 & 0 & \cdots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\0 & 0 & \cdots & 0 \end{bmatrix}. $$ The system of linear equations may be rewritten as: $$ (D+\omega L) \mathbf{x} = \omega \mathbf{b} - [\omega U + (\omega-1) D ] \mathbf{x} $$ for a constant \( \omega > 1 \).

The method of successive over-relaxation is an iterative technique that solves the left hand side of this expression for \( x \), using previous value for \( x \) on the right hand side. Analytically, this may be written as: $$ \mathbf{x}^{(k+1)} = (D+\omega L)^{-1} \big(\omega \mathbf{b} - [\omega U + (\omega-1) D ] \mathbf{x}^{(k)}\big). $$ However, by taking advantage of the triangular form of \( (D+\omega L) \), the elements of \( x^{(k+1)} \) can be computed sequentially using forward substitution: $$ x^{(k+1)}_i = (1-\omega)x^{(k)}_i + \frac{\omega}{a_{ii}} \left(b_i - \sum_{j > i} a_{ij}x^{(k)}_j - \sum_{j < i} a_{ij}x^{(k+1)}_j \right),\quad i=1,2,\ldots,n. $$ The choice of relaxation factor is not necessarily easy, and depends upon the properties of the coefficient matrix. For symmetric, positive-definite matrices it can be proven that \( 0 < \omega < 2 \) will lead to convergence, but we are generally interested in faster convergence rather than just convergence.

Cubic spline interpolation is among one of the mostly used methods for interpolating between data points where the arguments are organized as ascending series. In the library program we supply such a function, based on the so-called cubic spline method to be described below.

A spline function consists of polynomial pieces defined on subintervals. The different subintervals are connected via various continuity relations.

Assume we have at our disposal \( n+1 \) points \( x_0, x_1, \dots x_n \) arranged so that \( x_0 < x_1 < x_2 < \dots x_{n-1} < x_n \) (such points are called knots). A spline function \( s \) of degree \( k \) with \( n+1 \) knots is defined as follows

• On every subinterval \( [x_{i-1},x_i) \) s is a polynomial of degree \( \le k \).
• \( s \) has \( k-1 \) continuous derivatives in the whole interval \( [x_0,x_n] \).

As an example, consider a spline function of degree \( k=1 \) defined as follows $$ s(x)=\left\{\begin{array}{cc} s_0(x)=a_0x+b_0 & x\in [x_0, x_1) \\ s_1(x)=a_1x+b_1 & x\in [x_1, x_2) \\ \dots & \dots \\ s_{n-1}(x)=a_{n-1}x+b_{n-1} & x\in [x_{n-1}, x_n] \end{array}\right. $$ In this case the polynomial consists of series of straight lines connected to each other at every endpoint. The number of continuous derivatives is then \( k-1=0 \), as expected when we deal with straight lines. Such a polynomial is quite easy to construct given \( n+1 \) points \( x_0, x_1, \dots x_n \) and their corresponding function values.

The most commonly used spline function is the one with \( k=3 \), the so-called cubic spline function. Assume that we have in adddition to the \( n+1 \) knots a series of functions values \( y_0=f(x_0), y_1=f(x_1), \dots y_n=f(x_n) \). By definition, the polynomials \( s_{i-1} \) and \( s_i \) are thence supposed to interpolate the same point \( i \), that is $$ s_{i-1}(x_i)= y_i = s_i(x_i), $$ with \( 1 \le i \le n-1 \). In total we have \( n \) polynomials of the type $$ s_i(x)=a_{i0}+a_{i1}x+a_{i2}x^2+a_{i2}x^3, $$ yielding \( 4n \) coefficients to determine.

Every subinterval provides in addition the \( 2n \) conditions $$ y_i = s(x_i), $$ and $$ s(x_{i+1})= y_{i+1}, $$ to be fulfilled. If we also assume that \( s' \) and \( s'' \) are continuous, then $$ s'_{i-1}(x_i)= s'_i(x_i), $$ yields \( n-1 \) conditions. Similarly, $$ s''_{i-1}(x_i)= s''_i(x_i), $$ results in additional \( n-1 \) conditions. In total we have \( 4n \) coefficients and \( 4n-2 \) equations to determine them, leaving us with \( 2 \) degrees of freedom to be determined.

Using the last equation we define two values for the second derivative, namely $$ s''_{i}(x_i)= f_i, $$ and $$ s''_{i}(x_{i+1})= f_{i+1}, $$ and setting up a straight line between \( f_i \) and \( f_{i+1} \) we have $$ s_i''(x) = \frac{f_i}{x_{i+1}-x_i}(x_{i+1}-x)+ \frac{f_{i+1}}{x_{i+1}-x_i}(x-x_i), $$ and integrating twice one obtains $$ s_i(x) = \frac{f_i}{6(x_{i+1}-x_i)}(x_{i+1}-x)^3+ \frac{f_{i+1}}{6(x_{i+1}-x_i)}(x-x_i)^3 +c(x-x_i)+d(x_{i+1}-x). $$

Using the conditions \( s_i(x_i)=y_i \) and \( s_i(x_{i+1})=y_{i+1} \) we can in turn determine the constants \( c \) and \( d \) resulting in $$ \begin{eqnarray} s_i(x) =&\frac{f_i}{6(x_{i+1}-x_i)}(x_{i+1}-x)^3+ \frac{f_{i+1}}{6(x_{i+1}-x_i)}(x-x_i)^3 \nonumber \\ +&(\frac{y_{i+1}}{x_{i+1}-x_i}-\frac{f_{i+1}(x_{i+1}-x_i)}{6}) (x-x_i)+ (\frac{y_{i}}{x_{i+1}-x_i}-\frac{f_{i}(x_{i+1}-x_i)}{6}) (x_{i+1}-x). \end{eqnarray} $$

How to determine the values of the second derivatives \( f_{i} \) and \( f_{i+1} \)? We use the continuity assumption of the first derivatives $$ s'_{i-1}(x_i)= s'_i(x_i), $$ and set \( x=x_i \). Defining \( h_i=x_{i+1}-x_i \) we obtain finally the following expression $$ h_{i-1}f_{i-1}+2(h_{i}+h_{i-1})f_i+h_if_{i+1}= \frac{6}{h_i}(y_{i+1}-y_i)-\frac{6}{h_{i-1}}(y_{i}-y_{i-1}), $$ and introducing the shorthands \( u_i=2(h_{i}+h_{i-1}) \), \( v_i=\frac{6}{h_i}(y_{i+1}-y_i)-\frac{6}{h_{i-1}}(y_{i}-y_{i-1}) \), we can reformulate the problem as a set of linear equations to be solved through e.g., Gaussian elemination

Gaussian elimination $$ \left[\begin{array}{cccccccc} u_1 & h_1 &0 &\dots & & & & \\ h_1 & u_2 & h_2 &0 &\dots & & & \\ 0 & h_2 & u_3 & h_3 &0 &\dots & & \\ \dots& & \dots &\dots &\dots &\dots &\dots & \\ &\dots & & &0 &h_{n-3} &u_{n-2} &h_{n-2} \\ & && & &0 &h_{n-2} &u_{n-1} \end{array}\right] \left[\begin{array}{c} f_1 \\ f_2 \\ f_3\\ \dots \\ f_{n-2} \\ f_{n-1} \end{array} \right] = \left[\begin{array}{c} v_1 \\ v_2 \\ v_3\\ \dots \\ v_{n-2}\\ v_{n-1} \end{array} \right]. $$ Note that this is a set of tridiagonal equations and can be solved through only \( O(n) \) operations.

The functions supplied in the program library are spline and splint. In order to use cubic spline interpolation you need first to call

spline(double x[], double y[], int n, double yp1,  double yp2, double y2[])

This function takes as input \( x[0,..,n - 1] \) and \( y[0,..,n - 1] \) containing a tabulation \( y_i = f(x_i) \) with \( x_0 < x_1 < .. < x_{n - 1} \) together with the first derivatives of \( f(x) \) at \( x_0 \) and \( x_{n-1} \), respectively. Then the function returns \( y2[0,..,n-1] \) which contains the second derivatives of \( f(x_i) \) at each point \( x_i \). \( n \) is the number of points. This function provides the cubic spline interpolation for all subintervals and is called only once.

Thereafter, if you wish to make various interpolations, you need to call the function

splint(double x[], double y[], double y2a[], int n, double x, double *y)

which takes as input the tabulated values \( x[0,..,n - 1] \) and \( y[0,..,n - 1] \) and the output y2a[0,..,n - 1] from spline. It returns the value \( y \) corresponding to the point \( x \).

• Ordinary differential equations, Runge-Kutta method,chapter 8
• Ordinary differential equations with boundary conditions: one-variable equations to be solved by shooting and Green's function methods, chapter 9
• We can solve such equations by a finite difference scheme as well, turning the equation into an eigenvalue problem. Still one variable. Done in projects 1 and 2.
• If we have more than one variable, we need to solve partial differential equations, see Chapter 10

The material on differential equations is covered by chapters 8, 9 and 10. Project 3 deals with ordinary differential equations (the solar system). Reading assignment for week 39 is chapter 8 of lecture notes.

The order of the ODE refers to the order of the derivative on the left-hand side in the equation $$ \begin{equation} \frac{dy}{dt}=f(t,y). \end{equation} $$ This equation is of first order and \( f \) is an arbitrary function. A second-order equation goes typically like $$ \begin{equation} \frac{d^2y}{dt^2}=f(t,\frac{dy}{dt},y). \end{equation} $$ A well-known second-order equation is Newton's second law $$ \begin{equation} m\frac{d^2x}{dt^2}=-kx, \label{eq:newton} \end{equation} $$ where \( k \) is the force constant. ODE depend only on one variable

partial differential equations like the time-dependent Schr\"odinger equation $$ \begin{equation} i\hbar\frac{\partial \psi({\bf x},t)}{\partial t}= -\frac{\hbar^2}{2m}\left( \frac{\partial^2 \psi({\bf r},t)}{\partial x^2} + \frac{\partial^2 \psi({\bf r},t)}{\partial y^2}+ \frac{\partial^2 \psi({\bf r},t)}{\partial z^2}\right) + V({\bf x})\psi({\bf x},t), \end{equation} $$ may depend on several variables. In certain cases, like the above equation, the wave function can be factorized in functions of the separate variables, so that the Schroedinger equation can be rewritten in terms of sets of ordinary differential equations. These equations are discussed in chapter 10. Involve boundary conditions in addition to initial conditions.

We distinguish also between linear and non-linear differential equation where for example $$ \begin{equation} \frac{dy}{dt}=g^3(t)y(t), \end{equation} $$ is an example of a linear equation, while $$ \begin{equation} \frac{dy}{dt}=g^3(t)y(t)-g(t)y^2(t), \end{equation} $$ is a non-linear ODE.

Another concept which dictates the numerical method chosen for solving an ODE, is that of initial and boundary conditions. To give an example, if we study white dwarf stars or neutron stars we will need to solve two coupled first-order differential equations, one for the total mass \( m \) and one for the pressure \( P \) as functions of \( \rho \) $$ \frac{dm}{dr}=4\pi r^{2}\rho (r)/c^2, $$ and $$ \frac{dP}{dr}=-\frac{Gm(r)}{r^{2}}\rho (r)/c^2. $$ where \( \rho \) is the mass-energy density. The initial conditions are dictated by the mass being zero at the center of the star, i.e., when \( r=0 \), yielding \( m(r=0)=0 \). The other condition is that the pressure vanishes at the surface of the star.

In the solution of the Schroedinger equation for a particle in a potential, we may need to apply boundary conditions as well, such as demanding continuity of the wave function and its derivative.

In many cases it is possible to rewrite a second-order differential equation in terms of two first-order differential equations. Consider again the case of Newton's second law in Eq. \eqref{eq:newton}. If we define the position \( x(t)=y^{(1)}(t) \) and the velocity \( v(t)=y^{(2)}(t) \) as its derivative $$ \begin{equation} \frac{dy^{(1)}(t)}{dt}=\frac{dx(t)}{dt}=y^{(2)}(t), \end{equation} $$ we can rewrite Newton's second law as two coupled first-order differential equations $$ \begin{equation} m\frac{dy^{(2)}(t)}{dt}=-kx(t)=-ky^{(1)}(t), \label{eq:n1} \end{equation} $$ and $$ \begin{equation} \frac{dy^{(1)}(t)}{dt}=y^{(2)}(t). \label{eq:n2} \end{equation} $$

These methods fall under the general class of one-step methods. The algoritm is rather simple. Suppose we have an initial value for the function \( y(t) \) given by $$ \begin{equation} y_0=y(t=t_0). \end{equation} $$ We are interested in solving a differential equation in a region in space \( [a,b] \). We define a step \( h \) by splitting the interval in \( N \) sub intervals, so that we have $$ \begin{equation} h=\frac{b-a}{N}. \end{equation} $$ With this step and the derivative of \( y \) we can construct the next value of the function \( y \) at $$ \begin{equation} y_1=y(t_1=t_0+h), \end{equation} $$ and so forth.

If the function is rather well-behaved in the domain \( [a,b] \), we can use a fixed step size. If not, adaptive steps may be needed. Here we concentrate on fixed-step methods only. Let us try to generalize the above procedure by writing the step \( y_{i+1} \) in terms of the previous step \( y_i \) $$ \begin{equation} y_{i+1}=y(t=t_i+h)=y(t_i) + h\Delta(t_i,y_i(t_i)) + O(h^{p+1}), \end{equation} $$ where \( O(h^{p+1}) \) represents the truncation error. To determine \( \Delta \), we Taylor expand our function \( y \) $$ \begin{equation} y_{i+1}=y(t=t_i+h)=y(t_i) + h(y'(t_i)+\dots +y^{(p)}(t_i)\frac{h^{p-1}}{p!}) + O(h^{p+1}), \label{eq:taylor} \end{equation} $$ where we will associate the derivatives in the parenthesis with $$ \begin{equation} \Delta(t_i,y_i(t_i))=(y'(t_i)+\dots +y^{(p)}(t_i)\frac{h^{p-1}}{p!}). \label{eq:delta} \end{equation} $$

We define $$ \begin{equation} y'(t_i)=f(t_i,y_i) \end{equation} $$ and if we truncate \( \Delta \) at the first derivative, we have $$ \begin{equation} y_{i+1}=y(t_i) + hf(t_i,y_i) + O(h^2), \label{eq:euler} \end{equation} $$ which when complemented with \( t_{i+1}=t_i+h \) forms the algorithm for the well-known Euler method. Note that at every step we make an approximation error of the order of \( O(h^2) \), however the total error is the sum over all steps \( N=(b-a)/h \), yielding thus a global error which goes like \( NO(h^2)\approx O(h) \).

To make Euler's method more precise we can obviously decrease \( h \) (increase \( N \)). However, if we are computing the derivative \( f \) numerically by for example the two-steps formula $$ f'_{2c}(x)= \frac{f(x+h)-f(x)}{h}+O(h), $$ we can enter into roundoff error problems when we subtract two almost equal numbers \( f(x+h)-f(x)\approx 0 \). Euler's method is not recommended for precision calculation, although it is handy to use in order to get a first view on how a solution may look like. As an example, consider Newton's equation rewritten in Eqs. \eqref{eq:n1} and \eqref{eq:n2}. We define \( y_0=y^{(1)}(t=0) \) an \( v_0=y^{(2)}(t=0) \). The first steps in Newton's equations are then $$ \begin{equation} y^{(1)}_1=y_0+hv_0+O(h^2) \end{equation} $$ and $$ \begin{equation} y^{(2)}_1=v_0-hy_0k/m+O(h^2). \end{equation} $$

The Euler method is asymmetric in time, since it uses information about the derivative at the beginning of the time interval. This means that we evaluate the position at \( y^{(1)}_1 \) using the velocity at \( y^{(2)}_0=v_0 \). A simple variation is to determine \( y^{(1)}_{n+1} \) using the velocity at \( y^{(2)}_{n+1} \), that is (in a slightly more generalized form) $$ \begin{equation} y^{(1)}_{n+1}=y^{(1)}_{n}+h y^{(2)}_{n+1}+O(h^2) \end{equation} $$ and $$ \begin{equation} y^{(2)}_{n+1}=y^{(2)}_{n}+h a_{n}+O(h^2). \end{equation} $$ The acceleration \( a_n \) is a function of \( a_n(y^{(1)}_{n}, y^{(2)}_{n},t) \) and needs to be evaluated as well. This is the Euler-Cromer method.

Let us then include the second derivative in our Taylor expansion. We have then $$ \begin{equation} \Delta(t_i,y_i(t_i))=f(t_i)+\frac{h}{2}\frac{df(t_i,y_i)}{dt}+O(h^3). \end{equation} $$ The second derivative can be rewritten as $$ \begin{equation} y''=f'=\frac{df}{dt}=\frac{\partial f}{\partial t}+\frac{\partial f}{\partial y}\frac{\partial y}{\partial t}=\frac{\partial f}{\partial t}+\frac{\partial f}{\partial y}f \label{eq:derivatives} \end{equation} $$ and we can rewrite Eq.\ \eqref{eq:taylor} as $$ \begin{equation} y_{i+1}=y(t=t_i+h)=y(t_i) +hf(t_i)+ \frac{h^2}{2}\left(\frac{\partial f}{\partial t}+\frac{\partial f}{\partial y}f\right) + O(h^{3 }), \end{equation} $$ which has a local approximation error \( O(h^{3 }) \) and a global error \( O(h^{2}) \).

These approximations can be generalized by using the derivative \( f \) to arbitrary order so that we have $$ \begin{equation} y_{i+1}=y(t=t_i+h)=y(t_i) + h(f(t_i,y_i)+\dots f^{(p-1)}(t_i,y_i) \frac{h^{p-1}}{p!}) + O(h^{p+1}). \end{equation} $$ These methods, based on higher-order derivatives, are in general not used in numerical computation, since they rely on evaluating derivatives several times. Unless one has analytical expressions for these, the risk of roundoff errors is large.

The most obvious improvements to Euler's and Euler-Cromer's algorithms, avoiding in addition the need for computing a second derivative, is the so-called midpoint method. We have then $$ \begin{equation} y^{(1)}_{n+1}=y^{(1)}_{n}+\frac{h}{2}\left(y^{(2)}_{n+1}+y^{(2)}_{n}\right)+O(h^2) \end{equation} $$ and $$ \begin{equation} y^{(2)}_{n+1}=y^{(2)}_{n}+h a_{n}+O(h^2), \end{equation} $$ yielding $$ \begin{equation} y^{(1)}_{n+1}=y^{(1)}_{n}+hy^{(2)}_{n}+\frac{h^2}{2}a_n+O(h^3) \end{equation} $$ implying that the local truncation error in the position is now \( O(h^3) \), whereas Euler's or Euler-Cromer's methods have a local error of \( O(h^2) \).

Thus, the midpoint method yields a global error with second-order accuracy for the position and first-order accuracy for the velocity. However, although these methods yield exact results for constant accelerations, the error increases in general with each time step.

One method that avoids this is the so-called half-step method. Here we define $$ \begin{equation} y^{(2)}_{n+1/2}=y^{(2)}_{n-1/2}+h a_{n}+O(h^2), \end{equation} $$ and $$ \begin{equation} y^{(1)}_{n+1}=y^{(1)}_{n}+hy^{(2)}_{n+1/2} +O(h^2). \end{equation} $$ Note that this method needs the calculation of \( y^{(2)}_{1/2} \). This is done using e.g., Euler's method $$ \begin{equation} y^{(2)}_{1/2}=y^{(2)}_{0}+h a_{0}+O(h^2). \end{equation} $$ As this method is numerically stable, it is often used instead of Euler's method.

Another method which one may encounter is the Euler-Richardson method with $$ \begin{equation} y^{(2)}_{n+1}=y^{(2)}_{n}+h a_{n+1/2}+O(h^2), \label{eq:er1} \end{equation} $$ and $$ \begin{equation} \label{eq:er2} y^{(1)}_{n+1}=y^{(1)}_{n}+hy^{(2)}_{n+1/2} +O(h^2). \end{equation} $$ The program program2.cpp includes all of the above methods.

Runge-Kutta (RK) methods are based on Taylor expansion formulae, but yield in general better algorithms for solutions of an ODE. The basic philosophy is that it provides an intermediate step in the computation of \( y_{i+1} \).

To see this, consider first the following definitions $$ \begin{equation} \frac{dy}{dt}=f(t,y), \end{equation} $$ and $$ \begin{equation} y(t)=\int f(t,y) dt, \end{equation} $$ and $$ \begin{equation} y_{i+1}=y_i+ \int_{t_i}^{t_{i+1}} f(t,y) dt. \end{equation} $$

To demonstrate the philosophy behind RK methods, let us consider the second-order RK method, RK2. The first approximation consists in Taylor expanding \( f(t,y) \) around the center of the integration interval \( t_i \) to \( t_{i+1} \), that is, at \( t_i+h/2 \), \( h \) being the step. Using the midpoint formula for an integral, defining \( y(t_i+h/2) = y_{i+1/2} \) and \( t_i+h/2 = t_{i+1/2} \), we obtain $$ \begin{equation} \int_{t_i}^{t_{i+1}} f(t,y) dt \approx hf(t_{i+1/2},y_{i+1/2}) +O(h^3). \end{equation} $$ This means in turn that we have $$ \begin{equation} y_{i+1}=y_i + hf(t_{i+1/2},y_{i+1/2}) +O(h^3). \end{equation} $$

However, we do not know the value of \( y_{i+1/2} \). Here comes thus the next approximation, namely, we use Euler's method to approximate \( y_{i+1/2} \). We have then $$ \begin{equation} y_{(i+1/2)}=y_i + \frac{h}{2}\frac{dy}{dt} = y(t_i) + \frac{h}{2}f(t_i,y_i). \end{equation} $$ This means that we can define the following algorithm for the second-order Runge-Kutta method, RK2. $$ \begin{equation} k_1=hf(t_i,y_i), \end{equation} $$ $$ \begin{equation} k_2=hf(t_{i+1/2},y_i+k_1/2), \end{equation} $$ with the final value $$ \begin{equation} y_{i+i}\approx y_i + k_2 +O(h^3). \end{equation} $$

The difference between the previous one-step methods is that we now need an intermediate step in our evaluation, namely \( t_i+h/2 = t_{(i+1/2)} \) where we evaluate the derivative \( f \). This involves more operations, but the gain is a better stability in the solution.

The fourth-order Runge-Kutta, RK4, which we will employ in the solution of various differential equations below, has the following algorithm $$ \begin{equation} k_1=hf(t_i,y_i), \end{equation} $$ $$ \begin{equation} k_2=hf(t_i+h/2,y_i+k_1/2), \end{equation} $$ $$ \begin{equation} k_3=hf(t_i+h/2,y_i+k_2/2) \end{equation} $$ $$ \begin{equation} k_4=hf(t_i+h,y_i+k_3) \end{equation} $$ with the final value $$ \begin{equation} y_{i+1}=y_i +\frac{1}{6}\left( k_1 +2k_2+2k_3+k_4\right). \end{equation} $$ Thus, the algorithm consists in first calculating \( k_1 \) with \( t_i \), \( y_1 \) and \( f \) as inputs. Thereafter, we increase the step size by \( h/2 \) and calculate \( k_2 \), then \( k_3 \) and finally \( k_4 \). Global error as \( O(h^4) \).

Our first example is the classical case of simple harmonic oscillations, namely a block sliding on a horizontal frictionless surface. The block is tied to a wall with a spring. If the spring is not compressed or stretched too far, the force on the block at a given position \( x \) is $$ F=-kx. $$ The negative sign means that the force acts to restore the object to an equilibrium position. Newton's equation of motion for this idealized system is then $$ m\frac{d^2x}{dt^2}=-kx, $$ or we could rephrase it as $$ \frac{d^2x}{dt^2}=-\frac{k}{m}x=-\omega_0^2x, \label{eq:newton1} $$ with the angular frequency \( \omega_0^2=k/m \).

The above differential equation has the advantage that it can be solved analytically with solutions on the form $$ x(t)=Acos(\omega_0t+\nu), $$ where \( A \) is the amplitude and \( \nu \) the phase constant. This provides in turn an important test for the numerical solution and the development of a program for more complicated cases which cannot be solved analytically.

With the position \( x(t) \) and the velocity \( v(t)=dx/dt \) we can reformulate Newton's equation in the following way $$ \frac{dx(t)}{dt}=v(t), $$ and $$ \frac{dv(t)}{dt}=-\omega_0^2x(t). $$ We are now going to solve these equations using the Runge-Kutta method to fourth order discussed previously.

Before proceeding however, it is important to note that in addition to the exact solution, we have at least two further tests which can be used to check our solution.

Since functions like \( cos \) are periodic with a period \( 2\pi \), then the solution \( x(t) \) has also to be periodic. This means that $$ x(t+T)=x(t), $$ with \( T \) the period defined as $$ T=\frac{2\pi}{\omega_0}=\frac{2\pi}{\sqrt{k/m}}. $$ Observe that \( T \) depends only on \( k/m \) and not on the amplitude of the solution.

In addition to the periodicity test, the total energy has also to be conserved.

Suppose we choose the initial conditions $$ x(t=0)=1\hspace{0.1cm} \mathrm{m}\hspace{1cm} v(t=0)=0\hspace{0.1cm}\mathrm{m/s}, $$ meaning that block is at rest at \( t=0 \) but with a potential energy $$ E_0=\frac{1}{2}kx(t=0)^2=\frac{1}{2}k. $$ The total energy at any time \( t \) has however to be conserved, meaning that our solution has to fulfil the condition $$ E_0=\frac{1}{2}kx(t)^2+\frac{1}{2}mv(t)^2. $$

An algorithm which implements these equations is included below.

• Choose the initial position and speed, with the most common choice \( v(t=0)=0 \) and some fixed value for the position.
• Choose the method you wish to employ in solving the problem.
• Subdivide the time interval \( [t_i,t_f] \) into a grid with step size

$$ h=\frac{t_f-t_i}{N}, $$ where \( N \) is the number of mesh points.

• Calculate now the total energy given by

$$ E_0=\frac{1}{2}kx(t=0)^2=\frac{1}{2}k. $$

• The Runge-Kutta method is used to obtain \( x_{i+1} \) and \( v_{i+1} \) starting from the previous values \( x_i \) and \( v_i \).
• When we have computed \( x(v)_{i+1} \) we upgrade \( t_{i+1}=t_i+h \).
• This iterative process continues till we reach the maximum time \( t_f \).
• The results are checked against the exact solution. Furthermore, one has to check the stability of the numerical solution against the chosen number of mesh points \( N \).

/*    This program solves Newton's equation for a block
      sliding on a horizontal frictionless surface. The block
      is tied  to a wall with a spring, and Newton's equation
      takes the form
           m d^2x/dt^2 =-kx
      with k the spring tension and m the mass of the block.
      The angular frequency is omega^2 = k/m and we set it equal
      1 in this example program. 

      Newton's equation is rewritten as two coupled differential
      equations, one for the position x  and one for the velocity v
           dx/dt = v    and
           dv/dt = -x   when we set k/m=1

      We use therefore a two-dimensional array to represent x and v
      as functions of t
      y[0] == x
      y[1] == v
      dy[0]/dt = v
      dy[1]/dt = -x

      The derivatives are calculated by the user defined function 
      derivatives.

      The user has to specify the initial velocity (usually v_0=0)
      the number of steps and the initial position. In the programme
      below we fix the time interval [a,b] to [0,2*pi].

*/ 
#include <cmath>
#include <iostream>
#include <fstream>
#include <iomanip>
//#include "lib.h"
using namespace  std;
// output file as global variable
ofstream ofile;
// function declarations
void derivatives(double, double *, double *);
void initialise ( double&, double&, int&);
void output( double, double *, double);
void runge_kutta_4(double *, double *, int, double, double, 
                   double *, void (*)(double, double *, double *));

int main(int argc, char* argv[])
{
//  declarations of variables
  double *y, *dydt, *yout, t, h, tmax, E0;
  double initial_x, initial_v;
  int i, number_of_steps, n;
  char *outfilename;
  // Read in output file, abort if there are too few command-line arguments
  if( argc <= 1 ){
    cout << "Bad Usage: " << argv[0] <<
      " read also output file on same line" << endl;
    //    exit(1);
  }
  else{
    outfilename=argv[1];
  }
  ofile.open(outfilename);
  //  this is the number of differential equations  
  n = 2;     
  //  allocate space in memory for the arrays containing the derivatives 
  dydt = new double[n];
  y = new double[n];
  yout = new double[n];
  // read in the initial position, velocity and number of steps 
  initialise (initial_x, initial_v, number_of_steps);
  //  setting initial values, step size and max time tmax  
  h = 4.*acos(-1.)/( (double) number_of_steps);   // the step size     
  tmax = h*number_of_steps;               // the final time    
  y[0] = initial_x;                       // initial position  
  y[1] = initial_v;                       // initial velocity  
  t=0.;                                   // initial time      
  E0 = 0.5*y[0]*y[0]+0.5*y[1]*y[1];       // the initial total energy
  // now we start solving the differential equations using the RK4 method 
  while (t <= tmax){
    derivatives(t, y, dydt);   // initial derivatives              
    runge_kutta_4(y, dydt, n, t, h, yout, derivatives); 
    for (i = 0; i < n; i++) {
      y[i] = yout[i];  
    }
    t += h;
    output(t, y, E0);   // write to file 
  }
  delete [] y; delete [] dydt; delete [] yout;
  ofile.close();  // close output file
  return 0;
}   //  End of main function 

//     Read in from screen the number of steps,
//     initial position and initial speed 
void initialise (double& initial_x, double& initial_v, int& number_of_steps)
{
 cout << "Initial position = ";
 cin >> initial_x;
 cout << "Initial speed = ";
 cin >> initial_v;
 cout << "Number of steps = ";
 cin >> number_of_steps;
}  // end of function initialise  

//   this function sets up the derivatives for this special case  
void derivatives(double t, double *y, double *dydt)
{
  dydt[0]=y[1];    // derivative of x 
  dydt[1]=-y[0]; // derivative of v 
} // end of function derivatives  

//    function to write out the final results
void output(double t, double *y, double E0)
{
  ofile << setiosflags(ios::showpoint | ios::uppercase);
  ofile << setw(15) << setprecision(8) << t;
  ofile << setw(15) << setprecision(8) << y[0];
  ofile << setw(15) << setprecision(8) << y[1];
  ofile << setw(15) << setprecision(8) << cos(t);
  ofile << setw(15) << setprecision(8) << 
    0.5*y[0]*y[0]+0.5*y[1]*y[1]-E0 << endl;
}  // end of function output

/*   This function upgrades a function y (input as a pointer)
     and returns the result yout, also as a pointer. Note that
     these variables are declared as arrays.  It also receives as
     input the starting value for the derivatives in the pointer
     dydx. It receives also the variable n which represents the 
     number of differential equations, the step size h and 
     the initial value of x. It receives also the name of the
     function *derivs where the given derivative is computed
*/
void runge_kutta_4(double *y, double *dydx, int n, double x, double h, 
                   double *yout, void (*derivs)(double, double *, double *))
{
  int i;
  double      xh,hh,h6; 
  double *dym, *dyt, *yt;
  //   allocate space for local vectors   
  dym = new double [n];
  dyt =  new double [n];
  yt =  new double [n];
  hh = h*0.5;
  h6 = h/6.;
  xh = x+hh;
  for (i = 0; i < n; i++) {
    yt[i] = y[i]+hh*dydx[i];
  }
  (*derivs)(xh,yt,dyt);     // computation of k2, eq. 3.60   
  for (i = 0; i < n; i++) {
    yt[i] = y[i]+hh*dyt[i];
  }
  (*derivs)(xh,yt,dym); //  computation of k3, eq. 3.61   
  for (i=0; i < n; i++) {
    yt[i] = y[i]+h*dym[i];
    dym[i] += dyt[i];
  }
  (*derivs)(x+h,yt,dyt);    // computation of k4, eq. 3.62   
  //      now we upgrade y in the array yout  
  for (i = 0; i < n; i++){
    yout[i] = y[i]+h6*(dydx[i]+dyt[i]+2.0*dym[i]);
  }
  delete []dym;
  delete [] dyt;
  delete [] yt;
}       //  end of function Runge-kutta 4  

The angular equation of motion of the pendulum is given by Newton's equation and with no external force it reads $$ \begin{equation} ml\frac{d^2\theta}{dt^2}+mgsin(\theta)=0, \end{equation} $$ with an angular velocity and acceleration given by $$ \begin{equation} v=l\frac{d\theta}{dt}, \end{equation} $$ and $$ \begin{equation} a=l\frac{d^2\theta}{dt^2}. \end{equation} $$

We do however expect that the motion will gradually come to an end due a viscous drag torque acting on the pendulum. In the presence of the drag, the above equation becomes $$ \begin{equation} ml\frac{d^2\theta}{dt^2}+\nu\frac{d\theta}{dt} +mgsin(\theta)=0, \label{eq:pend1} \end{equation} $$ where \( \nu \) is now a positive constant parameterizing the viscosity of the medium in question. In order to maintain the motion against viscosity, it is necessary to add some external driving force. We choose here a periodic driving force. The last equation becomes then $$ \begin{equation} ml\frac{d^2\theta}{dt^2}+\nu\frac{d\theta}{dt} +mgsin(\theta)=Asin(\omega t), \label{eq:pend2} \end{equation} $$ with \( A \) and \( \omega \) two constants representing the amplitude and the angular frequency respectively. The latter is called the driving frequency.

We define $$ \omega_0=\sqrt{g/l}, $$ the so-called natural frequency and the new dimensionless quantities $$ \hat{t}=\omega_0t, $$ with the dimensionless driving frequency $$ \hat{\omega}=\frac{\omega}{\omega_0}, $$ and introducing the quantity \( Q \), called the quality factor, $$ Q=\frac{mg}{\omega_0\nu}, $$ and the dimensionless amplitude $$ \hat{A}=\frac{A}{mg} $$

We have $$ \frac{d^2\theta}{d\hat{t}^2}+\frac{1}{Q}\frac{d\theta}{d\hat{t}} +sin(\theta)=\hat{A}cos(\hat{\omega}\hat{t}). $$ This equation can in turn be recast in terms of two coupled first-order differential equations as follows $$ \frac{d\theta}{d\hat{t}}=\hat{v}, $$ and $$ \frac{d\hat{v}}{d\hat{t}}=-\frac{\hat{v}}{Q}-sin(\theta)+\hat{A}cos(\hat{\omega}\hat{t}). $$ These are the equations to be solved. The factor \( Q \) represents the number of oscillations of the undriven system that must occur before its energy is significantly reduced due to the viscous drag. The amplitude \( \hat{A} \) is measured in units of the maximum possible gravitational torque while \( \hat{\omega} \) is the angular frequency of the external torque measured in units of the pendulum's natural frequency.

It can be very useful to make a Class which contains all possible methods discussed. In Fortran we can use the MODULE keyword in order to can methods and keep the variables private and hidden from other parts of our code. This allows for a generalization which can be used to tackle other ODEs as well. In program2.cpp of chapter 8 we have canned the following methods

• void euler();
• void euler_cromer();
• void midpoint();
• void euler_richardson();
• void half_step();
• void rk2(); //runge-kutta-second-order
• void rk4_step(double,double*,double*,double); // we need it in function rk4() and asc()
• void rk4(); //runge-kutta-fourth-order
• void asc(); //runge-kutta-fourth-order with adaptive stepsize control

#include <stdio.h>
#include <iostream>
#include <cmath>
#include <fstream>
#include <iomanip>
using namespace  std;
/*

Different methods for solving ODEs are presented
We are solving the following eqation:

m*l*(phi)'' + reib*(phi)' + m*g*sin(phi) = A*cos(omega*t)

If you want to solve similar equations with other values you have to
rewrite the methods 'derivatives' and 'initialise' and change the variables in the private
part of the class Pendel

At first we rewrite the equation using the following definitions:

omega_0 = sqrt(g*l)
t_roof = omega_0*t
omega_roof = omega/omega_0
Q = (m*g)/(omega_0*reib)
A_roof = A/(m*g)

and we get a dimensionless equation

(phi)'' + 1/Q*(phi)' + sin(phi) = A_roof*cos(omega_roof*t_roof)

This equation can be written as two equations of first order:

(phi)' = v
(v)' = -v/Q - sin(phi) +A_roof*cos(omega_roof*t_roof)

All numerical methods are applied to the last two equations.
*/

class pendelum
 {
 private:
   double Q, A_roof, omega_0, omega_roof,g; //
   double y[2];          //for the initial-values of phi and v
   int n;                // how many steps
   double delta_t,delta_t_roof;

 public:
   void derivatives(double,double*,double*);
   void initialise();
   void euler();
   void euler_cromer();
   void midpoint();
   void euler_richardson();
   void half_step();
   void rk2(); //runge-kutta-second-order
   void rk4_step(double,double*,double*,double); // we need it in function rk4() and asc()
   void rk4(); //runge-kutta-fourth-order
   void asc(); //runge-kutta-fourth-order with adaptive stepsize control
 };

void pendelum::derivatives(double t, double* in, double* out)
{ /* Here we are calculating the derivatives at (dimensionless) time t
     'in' are the values of phi and v, which are used for the calculation
     The results are given to 'out' */
  
  out[0]=in[1];             //out[0] = (phi)'  = v
  if(Q)
    out[1]=-in[1]/((double)Q)-sin(in[0])+A_roof*cos(omega_roof*t);  //out[1] = (phi)''
  else
    out[1]=-sin(in[0])+A_roof*cos(omega_roof*t);  //out[1] = (phi)''
}


void pendelum::initialise()
{
  double m,l,omega,A,viscosity,phi_0,v_0,t_end;
  cout<<"Solving the differential eqation of the pendulum!\n";
  cout<<"We have a pendulum with mass m, length l. Then we have a periodic force with amplitude A and omega\n";
  cout<<"Furthermore there is a viscous drag coefficient.\n";
  cout<<"The initial conditions at t=0 are phi_0 and v_0\n";
  cout<<"Mass m: ";
  cin>>m;
  cout<<"length l: ";
  cin>>l;
  cout<<"omega of the force: ";
  cin>>omega;
  cout<<"amplitude of the force: ";
  cin>>A;
  cout<<"The value of the viscous drag constant (viscosity): ";
  cin>>viscosity;
  cout<<"phi_0: ";
  cin>>y[0];
  cout<<"v_0: ";
  cin>>y[1];
  cout<<"Number of time steps or integration steps:";
  cin>>n;
  cout<<"Final time steps as multiplum of pi:";
  cin>>t_end;
  t_end *= acos(-1.);
  g=9.81;
  // We need the following values:
  omega_0=sqrt(g/((double)l));      // omega of the pendulum
  if (viscosity)  Q= m*g/((double)omega_0*viscosity);
  else Q=0; //calculating Q
  A_roof=A/((double)m*g);
  omega_roof=omega/((double)omega_0);
  delta_t_roof=omega_0*t_end/((double)n);    //delta_t without dimension
  delta_t=t_end/((double)n);
}

void pendelum::euler()
{ //using simple euler-method
  int i;
  double yout[2],y_h[2];
  double t_h;

  y_h[0]=y[0];
  y_h[1]=y[1];
  t_h=0;
  ofstream fout("euler.out");
  fout.setf(ios::scientific);
  fout.precision(20);
  for(i=1;i<=n;i++){
    derivatives(t_h,y_h,yout);
    yout[1]=y_h[1]+yout[1]*delta_t_roof;
    yout[0]=y_h[0]+yout[0]*delta_t_roof;
    // Calculation with dimensionless values    
    fout<<i*delta_t<<"\t\t"<<yout[0]<<"\t\t"<<yout[1]<<"\n";
    t_h+=delta_t_roof;
    y_h[1]=yout[1];
    y_h[0]=yout[0];
  }
  fout.close;
}

void pendelum::euler_cromer()
{
  int i;
  double t_h;
  double yout[2],y_h[2];

  t_h=0;
  y_h[0]=y[0];  //phi
  y_h[1]=y[1];  //v
  ofstream fout("ec.out");
  fout.setf(ios::scientific);
  fout.precision(20);
  for(i=1; i<=n; i++){
    derivatives(t_h,y_h,yout);
    yout[1]=y_h[1]+yout[1]*delta_t_roof;
    yout[0]=y_h[0]+yout[1]*delta_t_roof;
    // The new calculated value of v is used for calculating phi 
    fout<<i*delta_t<<"\t\t"<<yout[0]<<"\t\t"<<yout[1]<<"\n";
    t_h+=delta_t_roof;
    y_h[0]=yout[0];
    y_h[1]=yout[1];
  }
  fout.close;
}

void pendelum::midpoint()
{
  int i;
  double t_h;
  double yout[2],y_h[2];
  
  t_h=0;
  y_h[0]=y[0];  //phi
  y_h[1]=y[1];  //v
  ofstream fout("midpoint.out");
  fout.setf(ios::scientific);
  fout.precision(20);
  for(i=1; i<=n; i++){
    derivatives(t_h,y_h,yout);
    yout[1]=y_h[1]+yout[1]*delta_t_roof;
    yout[0]=y_h[0]+0.5*(yout[1]+y_h[1])*delta_t_roof;
    fout<<i*delta_t<<"\t\t"<<yout[0]<<"\t\t"<<yout[1]<<"\n";
    t_h+=delta_t_roof;
    y_h[0]=yout[0];
    y_h[1]=yout[1];
  }
  fout.close;
}


void pendelum::euler_richardson()
{
  int i;
  double t_h,t_m;
  double yout[2],y_h[2],y_m[2];

  t_h=0;
  y_h[0]=y[0];  //phi
  y_h[1]=y[1];  //v
  ofstream fout("er.out");
  fout.setf(ios::scientific);
  fout.precision(20);
  for(i=1; i<=n; i++){
    derivatives(t_h,y_h,yout);
    y_m[1]=y_h[1]+0.5*yout[1]*delta_t_roof;
    y_m[0]=y_h[0]+0.5*y_h[1]*delta_t_roof;
    t_m=t_h+0.5*delta_t_roof;
    derivatives(t_m,y_m,yout);
    yout[1]=y_h[1]+yout[1]*delta_t_roof;
    yout[0]=y_h[0]+y_m[1]*delta_t_roof;
    fout<<i*delta_t<<"\t\t"<<yout[0]<<"\t\t"<<yout[1]<<"\n";
    t_h+=delta_t_roof;
    y_h[0]=yout[0];
    y_h[1]=yout[1];
  }
  fout.close;
}

void pendelum::half_step()
{
  /*We are using the half_step_algorith.
    The algorithm is not self-starting, so we calculate
    v_1/2 by using the Euler algorithm. */

  int i;
  double t_h;
  double yout[2],y_h[2];

  t_h=0;
  y_h[0]=y[0];  //phi
  y_h[1]=y[1];  //v
  ofstream fout("half_step.out");
  fout.setf(ios::scientific);
  fout.precision(20);
  /*At first we have to calculate v_1/2
    For this we use Euler's method:
    v_`1/2 = v_0 + 1/2*a_0*delta_t_roof
    For calculating a_0 we have to start derivatives
  */
  derivatives(t_h,y_h,yout);
  yout[1]=y_h[1]+0.5*yout[1]*delta_t_roof;
  yout[0]=y_h[0]+yout[1]*delta_t_roof;
  fout<<delta_t<<"\t\t"<<yout[0]<<"\t\t"<<yout[1]<<"\n";
  y_h[0]=yout[0];
  y_h[1]=yout[1];
  for(i=2; i<=n; i++){
    derivatives(t_h,y_h,yout);
    yout[1]=y_h[1]+yout[1]*delta_t_roof;
    yout[0]=y_h[0]+yout[1]*delta_t_roof;
    fout<<i*delta_t<<"\t\t"<<yout[0]<<"\t\t"<<yout[1]<<"\n";
    t_h+=delta_t_roof;
    y_h[0]=yout[0];
    y_h[1]=yout[1];
  }
  fout.close;
}

void pendelum::rk2()
{
  /*We are using the second-order-Runge-Kutta-algorithm
    We have to calculate the parameters k1 and k2 for v and phi,
    so we use to arrays k1[2] and k2[2] for this
    k1[0], k2[0] are the parameters for phi,
    k1[1], k2[1] are the parameters for v
  */

  int i;
  double t_h;
  double yout[2],y_h[2],k1[2],k2[2],y_k[2];
  
  t_h=0;
  y_h[0]=y[0];  //phi
  y_h[1]=y[1];  //v
  ofstream fout("rk2.out");
  fout.setf(ios::scientific);
  fout.precision(20);
  for(i=1; i<=n; i++){
    /*Calculation of k1 */
    derivatives(t_h,y_h,yout);
    k1[1]=yout[1]*delta_t_roof;
    k1[0]=yout[0]*delta_t_roof;
    y_k[0]=y_h[0]+k1[0]*0.5;
    y_k[1]=y_h[1]+k2[1]*0.5;
    /*Calculation of k2 */
    derivatives(t_h+delta_t_roof*0.5,y_k,yout);
    k2[1]=yout[1]*delta_t_roof;
    k2[0]=yout[0]*delta_t_roof;
    yout[1]=y_h[1]+k2[1];
    yout[0]=y_h[0]+k2[0];
    fout<<i*delta_t<<"\t\t"<<yout[0]<<"\t\t"<<yout[1]<<"\n";
    t_h+=delta_t_roof;
    y_h[0]=yout[0];
    y_h[1]=yout[1];
  }
  fout.close;
}

void pendelum::rk4_step(double t,double *yin,double *yout,double delta_t)
{
  /*
    The function calculates one step of fourth-order-runge-kutta-method
    We will need it for the normal fourth-order-Runge-Kutta-method and
    for RK-method with adaptive stepsize control

    The function calculates the value of y(t + delta_t) using fourth-order-RK-method
    Input: time t and the stepsize delta_t, yin (values of phi and v at time t)
    Output: yout (values of phi and v at time t+delta_t)
    
  */
  double k1[2],k2[2],k3[2],k4[2],y_k[2];
  // Calculation of k1 
  derivatives(t,yin,yout);
  k1[1]=yout[1]*delta_t;
  k1[0]=yout[0]*delta_t;
  y_k[0]=yin[0]+k1[0]*0.5;
  y_k[1]=yin[1]+k1[1]*0.5;
  /*Calculation of k2 */
  derivatives(t+delta_t*0.5,y_k,yout);
  k2[1]=yout[1]*delta_t;
  k2[0]=yout[0]*delta_t;
  y_k[0]=yin[0]+k2[0]*0.5;
  y_k[1]=yin[1]+k2[1]*0.5;
  /* Calculation of k3 */
  derivatives(t+delta_t*0.5,y_k,yout);
  k3[1]=yout[1]*delta_t;
  k3[0]=yout[0]*delta_t;
  y_k[0]=yin[0]+k3[0];
  y_k[1]=yin[1]+k3[1];
  /*Calculation of k4 */
  derivatives(t+delta_t,y_k,yout);
  k4[1]=yout[1]*delta_t;
  k4[0]=yout[0]*delta_t;
  /*Calculation of new values of phi and v */
  yout[0]=yin[0]+1.0/6.0*(k1[0]+2*k2[0]+2*k3[0]+k4[0]);
  yout[1]=yin[1]+1.0/6.0*(k1[1]+2*k2[1]+2*k3[1]+k4[1]);
}

void pendelum::rk4()
{
  /*We are using the fourth-order-Runge-Kutta-algorithm
    We have to calculate the parameters k1, k2, k3, k4 for v and phi,
    so we use to arrays k1[2] and k2[2] for this
    k1[0], k2[0] are the parameters for phi,
    k1[1], k2[1] are the parameters for v
  */

  int i;
  double t_h;
  double yout[2],y_h[2]; //k1[2],k2[2],k3[2],k4[2],y_k[2];
  
  t_h=0;
  y_h[0]=y[0];  //phi
  y_h[1]=y[1];  //v
  ofstream fout("rk4.out");
  fout.setf(ios::scientific);
  fout.precision(20);
  for(i=1; i<=n; i++){
    rk4_step(t_h,y_h,yout,delta_t_roof);
    fout<<i*delta_t<<"\t\t"<<yout[0]<<"\t\t"<<yout[1]<<"\n";
    t_h+=delta_t_roof;
    y_h[0]=yout[0];
    y_h[1]=yout[1];
  }
  fout.close;
}



void pendelum::asc()
{
  /*
    We are using the Runge-Kutta-algorithm with adaptive stepsize control
    according to "Numerical Recipes in C", S. 574 ff.
    
    At first we calculate y(x+h) using rk4-method  => y1
    Then we calculate y(x+h) using two times rk4-method at x+h/2 and x+h  => y2
    
    The difference between these values is called "delta" If it is smaller than a given value,
    we calculate y(x+h) by  y2 + (delta)/15 (page 575, Numerical R.)
    
    If delta is not smaller than ... we calculate a new stepsize using
    h_new=(Safety)*h_old*(.../delta)^(0.25) where "Safety" is constant (page 577 N.R.)
    and start again with calculating y(x+h)...
   */
  int i;
  double t_h,h_alt,h_neu,hh,errmax;
  double yout[2],y_h[2],y_m[2],y1[2],y2[2], delta[2], yscal[2];
  
  const double eps=1.0e-6;
  const double safety=0.9;
  const double errcon=6.0e-4;
  const double tiny=1.0e-30;

  t_h=0;
  y_h[0]=y[0];  //phi
  y_h[1]=y[1];  //v
  h_neu=delta_t_roof;
  ofstream fout("asc.out");
  fout.setf(ios::scientific);
  fout.precision(20);
  for(i=0;i<=n;i++){
    /* The error is scaled against yscal
       We use a yscal of the form yscal = fabs(y[i]) + fabs(h*derivatives[i])
       (N.R. page 567)
    */
    derivatives(t_h,y_h,yout);
    yscal[0]=fabs(y[0])+fabs(h_neu*yout[0])+tiny;
    yscal[1]=fabs(y[1])+fabs(h_neu*yout[1])+tiny;
    /* the do-while-loop is used until the */
    do{
      /* Calculating y2 by two half steps */
      h_alt=h_neu;
      hh=h_alt*0.5;
      rk4_step(t_h, y_h, y_m, hh);
      rk4_step(t_h+hh,y_m,y2,hh);
      /* Calculating y1 by one normal step */
      rk4_step(t_h,y_h,y1,h_alt);
      /* Now we have two values for phi and v at the time t_h + h  in y2 and y1
	 We can now calculate the delta for phi and v
      */
      delta[0]=fabs(y1[0]-y2[0]);
      delta[1]=fabs(y1[1]-y2[1]);
      errmax=(delta[0]/yscal[0] > delta[1]/yscal[1] ? delta[0]/yscal[0] : delta[1]/yscal[1]);
      
      /*We scale delta against the constant yscal
	Then we take the biggest one and call it errmax */
      errmax=(double)errmax/eps;
      /*We divide errmax by eps and have only   */
      h_neu=safety*h_alt*exp(-0.25*log(errmax));
    }while(errmax>1.0);
    /*Now we are outside the do-while-loop and have a delta which is small enough
      So we can calculate the new values of phi and v
    */
    yout[0]=y_h[0]+delta[0]/15.0;
    yout[1]=y_h[1]+delta[1]/15.0;
    fout<<(double)(t_h+h_alt)/omega_0<<"\t\t"<<yout[0]<<"\t\t"<<yout[1]<<"\n";
    // Calculating of the new stepsize
    h_neu=(errmax > errcon ? safety*h_alt*exp(-0.20*log(errmax)) : 4.0*h_alt);
    y_h[0]=yout[0];
    y_h[1]=yout[1];
    t_h+=h_neu;
  }
}


int main()
{
  pendelum testcase;
  testcase.initialise();
  //  testcase.euler();
  //testcase.euler_cromer();
  //testcase.midpoint();
  //testcase.euler_richardson();
  //testcase.half_step();
  //testcase.rk2();
  testcase.rk4();
  return 0;
}  // end of main function

In Fortran we would use

MODULE pendulum
   USE CONSTANTS
   IMPLICIT NONE
   REAL(DP), PRIVATE :: Q, A_roof, omega_0, omega_roof,g
   REAL(DP), PRIVATE :: y(2)         ! for the initial-values of phi and v
   INTEGER, PRIVATE ::  n               ! how many steps
   REAL(DP), PRIVATE :: delta_t,delta_t_roof

   CONTAINS
    SUBROUTINE derivatives(..)
    SUBROUTINE initialise(..)
    SUBROUTINE euler(..)
    SUBROUTINE euler_cromer(..)
    SUBROUTINE midpoint(..)
    etc

END MODULE pendulum

In case the function to integrate varies slowly or fast in different integration domains, adaptive methods are normally used. One strategy is always to decrease the step size. As we have seen earlier, this leads to more CPU cycles and may lead to loss or numerical precision. An alternative is to use higher-order RK methods for example. However, this leads again to more cycles, furthermore, there is no guarantee that higher-order leads to an improved error.

Assume the exact result is \( \tilde{x} \) and that we are using an RKM method. Suppose we run two calculations, one with \( h \) (called \( x_1 \)) and one with \( h/2 \) (called \( x_2 \)). Then $$ \tilde{x}=x_1+Ch^{M+1}+O(h^{M+2}), $$ and $$ \tilde{x}=x_2+2C(h/2)^{M+1}+O(h^{M+2}), $$ with \( C \) a constant. Note that we calculate two halves in the last equation. We get then $$ |x_1-x_2| = Ch^{M+1}(1-\frac{1}{2^M}). $$ yielding $$ C=\frac{|x_1-x_2|}{(1-2^{-M})h^{M+1}}. $$ We rewrite $$ \tilde{x}=x_2+\epsilon+O((h)^{M+2}), $$ with $$ \epsilon = \frac{|x_1-x_2|}{2^M-1}. $$

With RK4 the expressions become $$ \tilde{x}=x_2+\epsilon+O((h)^{6}), $$ with $$ \epsilon = \frac{|x_1-x_2|}{15}. $$ The estimate is one order higher than the original RK4. But this method is normally rather inefficient since it requires a lot of computations. We solve typically the equation three times at each time step. However, we can compare the estimate \( \epsilon \) with some by us given accuracy \( \xi \). We can then ask the question: what is, with a given \( x_j \) and \( t_j \), the largest possible step size \( \tilde{h} \) that leads to a truncation error below \( \xi \)? We want $$ C\tilde{h} \le \xi, $$ which leads to $$ \left(\frac{\tilde{h}}{h}\right)^{M+1}\frac{|x_1-x_2|}{(1-2^{-M})}\le \xi, $$ meaning that $$ \tilde{h}=h\left(\frac{\xi}{\epsilon}\right)^{1+1/M}. $$

With $$ \tilde{h}=h\left(\frac{\xi}{\epsilon}\right)^{1+1/M}. $$ we can design the following algorithm:

• If the two answers are close, keep the approximation to \( h \).
• If \( \epsilon > \xi \) we need to decrease the step size in the next time step.
• If \( \epsilon < \xi \) we need to increase the step size in the next time step.

A much used algorithm is the so-called RKF45 which uses a combination of fourth and fifth order RK methods.

At each step, two different approximations for the solution are made and compared. If the two answers are in close agreement, the approximation is accepted. If the two answers do not agree to a specified accuracy, the step size is reduced. If the answers agree to more significant digits than required, the step size is increased. Each step requires the use of the following six values: $$ k_1 = h f (t_k , y_k ), $$ $$ k_2 = h f (t_k + \frac{1}{4}h, y_k + \frac{1}{4}k_1) , $$ $$ k_3 = h f (t_k + \frac{3}{8}h, y_k + \frac{3}{32}k_1 + \frac{9}{32}k_2) , $$ $$ k_4 = h f (t_k + \frac{12}{13}h, y_k + \frac{1932}{2197}k_1 + \frac{7200}{2197}k_2+\frac{7296}{2197}k_3), $$ $$ k_5 = h f (t_k + h, y_k + \frac{439}{216}k_1 -8k_2+ \frac{3680}{513}k_3+\frac{845}{4104}k_4), $$ $$ k_6 = h f (t_k + \frac{1}{2}h, y_k - \frac{8}{27}k_1 + 2k_2-\frac{3544}{2565}k_2+\frac{1859}{4104}k_4-+\frac{11}{40}k_5). $$

An approximation to the solution of the ODE is made using a Runge-Kutta method of order 4: $$ y_{k+1} = y_k + \frac{25}{216}k_1+\frac{1408}{2565}k_3 +\frac{2197}{4101}k_4-\frac{1}{5}k_5, $$ where the four function values \( k_1 \) , \( k_3 \) , \( k_4 \) , and \( k_5 \) are used. Notice that \( k_2 \) is not used here. A better value for the solution is determined using a Runge-Kutta method of order 5: $$ z_{k+1} = y_k + \frac{16}{135}k_1+\frac{6656}{12825}k_3 +\frac{28561}{56430}k_4-\frac{9}{50}k_5+\frac{2}{55}k_6. $$ The optimal time step \( \alpha h \) is then determined by $$ \alpha = \left( \frac{\xi h}{2|z_{k+1}-y_{k+1}|}\right)^{1/4}, $$ with \( \xi \) our defined tolerance.

See https://github.com/htihle/Keplerproblem. To be discussed during the lab sessions in week 40